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We use Monte Carlo simulations to study properties of Anderson's resonating-valence-bond (RVB) 
spin- liquid state on the square lattice (i.e., the equal superposition of all pairing of spins into nearest- 
neighbor singlet pairs) and compare with the classical dimer model (CDM). The latter system also 
corresponds to the ground state of the Rokhsar-Kivelson quantum dimer model at its critical point. 
We find that although spin-spin correlations decay exponentially in the RVB, four-spin valence- 
bond-solid correlations are critical, qualitatively like the well-known dimer-dimer correlations of the 
CDM, but decaying more slowly (as l/r a with a m 1.20, compared with a — 2 for the CDM). 
We also compute the distribution of monomer (defect) pair separations, which decay by a larger 
exponent in the RVB than in the CDM. We further study both models in their different winding 
number sectors and evaluate the relative weights of different sectors. Like the CDM, all the observed 
RVB behaviors can be understood in the framework of a mapping to a "height" model characterized 
by a gradient-squared stiffness constant K. Four independent measurements consistently show a 
value A'rvb ~ 1.6-Kcdm, with the same kinds of numerical evaluations of A'cdm giving results in 
agreement with the rigorously known value A'cdm = 7r/16. The background of a nonzero winding 
number gradient W/L introduces spatial anisotropics and an increase in the effective K, both of 
which can be understood as a consequence of anharmonic terms in the height-model free energy, 
which are of relevance to the recently proposed scenario of "Cantor deconfinement" in extended 
quantum dimer models. In addition to the standard case of short bonds only, we also studied 
ensembles in which fourth-neighbor (bipartite) bonds are allowed, at a density controlled by a 
tunable fugacity, resulting (as expected) in a smooth reduction of K. 

PACS numbers: 75.10.Jm, 75.10.Nr, 75.40.Mg, 75.40.Cx 



I. INTRODUCTION 

The two-dimensional (2D) resonating-valence-bond 
(RVB) spin-liquid state introduced by Anderson has been 
studied extensively during the past two decades, with 
the hope that it (when doped) might provide an oppor- 
tunity to understand high-temperature superconductiv- 
ity in cupratesi Such RVB states, which do not feature 
any long range magnetic order or broken lattice symme- 
tries (but are believed to exhibit non-local, topological 
order— are also of broader interest in the context of 
frustrated magnetism, where they were first considered^ 
In studies of specific Hamiltonians, RVB states can be 
considered as variational ground states. The extreme 
RVB state built out of only the shortest possible (nearest- 
neighbor) valence bonds (singlets), with equal weights 
for all bond configurations (which in the case considered 
here will be on the square lattice), does not have any 
adjustable parameters (as long as the signs of the wave 
function are not considered — in the standard RVB all 
coefficients are equal and positive). One can also para- 
metrically introduce longer bonds in amplitude-product 
states^ In two dimensions these states are spin liquids if 
the amplitudes decay sufficiently rapidly (exponentially 
or as a high power) with the bond length. We report 
here extensive studies of the RVB state, with only short 
(length 1) bonds, as well as in the presence of a fraction 
of bonds (the second bipartite ones of length . 



The search for Hamiltonians with RVB ground states 
has been an ongoing challenge during the past two 
decades. One way to approach the problem is through 
quantum dimer models (QDM), in which the internal sin- 
glet structure of the valence bonds is neglected. The 
valence bonds are replaced by hard-core dimers, and dif- 
ferent dimer configurations are considered as orthogonal 
states^ The effective Hamiltonians in this space, which 
describe the quantum fluctuations of the dimers, can have 
crystalline dimer order [corresponding to a valcncc-bond- 
solid (VBS) in the spin system] or be disordered (corre- 
sponding to a spin liquid). QDMs have many interest- 
ing and intriguing properties, e.g., the special Rokhsar- 
Kivelson (RK) points at which the wave-function of a 
dimer model corresponds exactly to the statistical me- 
chanics of classical dimers.— ~— On the square lattice the 
classical dimer model (CDM) has critical dimer-dimer 
correlations, decaying with distance r as 1/r 2 (a rig- 
orous result^) which then is also the case at the RK 
point separating two different VBS states on the square 
lattice. On the triangular lattice, this isolated spin- 
liquid point with critical dimer correlations is replaced 
by an extended liquid phase with exponentially decaying 
dimer correlations ^ The same physics can be achieved 
on the square lattice by introducing dimers between next- 
nearest-neighbor sitesJ^ We will here also provide some 
further results for the CDM, in order to elucidate in more 
detail the relationship between the RVB and the CDM. 
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Formally, the QDMs can be related exactly to gen- 
eralized SU(iV) symmetric spin models^ In the limit 
of N — s- oo the valence-bond states become exactly or- 
thogonal dimer states. Whether or not the physics of 
the quantum dimer models can be extended down to the 
physically most interesting case of SU(2) spins is in gen- 
eral not clear (unless the N = 2 features are built in 
from the start, as can be done in generalized QDMs^-). 
Moessner and Sondhi have devised a procedure to mimic 
a system of large- N spins by decorating an original lat- 
tice of S = 1/2 SU(2) spins with additional spins, and 
this way a Hamiltonian with spin-liquid ground state 
can be constructed^ Very recently, Cano and Fendley 
constructed a Hamiltonian the ground state of which is 
exactly the short-bond RVB state on the square lattice 
(without decoration)^ While this Hamiltonian is a com- 
plicated one with multi-spin interactions that are unlikely 
present in real systems, the achievement is important as 
it shows that local SU(2) spin models with RVB states 
do in principle exist also on simple lattices. 

A. Correlations in RVB and dimer states 

Perhaps surprisingly, very few physical properties of 
RVB spin liquids have actually been computed. While 
Monte Carlo simulations of amplitude-product states on 
the 2D square lattice were carried out some time ago, 
only the simple spin-spin correlations were calculated.— 
They decay exponentially in the case of the short-bond 
state. On the other hand, the fact that the dimer-dimcr 
correlations of the CDM (or, cquivalcntly, the QDM at 
the RK point) decay with a power-law clearly suggests 
that there should be similar critical correlations also in 
the RVB state (if the QDM is qualitatively faithful to 
it). The dimer-dimer correlations of the RVB state are 
not physical correlations, however, as the dimer basis is 
non-orthogonal and ovcrcompletc. 

In this paper, we use an improved Monte Carlo sam- 
pling scheme for valence bondsii to compute the physical 
correlation function most closely related to the dimer- 
dimer correlations of the CDM, namely, the four-spin 
correlation function 

Acxfc,) = (BxfojBxfa)), (1.1) 

where B x (y{) is a scalar operator defined on a bond, 

B x (n) = S(n) ■ S(n + ft), (1.2) 

and D yy and D xy can be defined analogously. Here the 
lattice coordinate of spin i is denoted and x is the 
lattice vector in the x-direction. The operator B x (r.i) 
provides a measure of the singlet probability on the bond 
between site i and its "right" neighbor, which is larger on 
a valence bond (in which case the operator is diagonal) 
than between two valence bonds (where the operator is 
off-diagonal and leads to a rearrangement of the two va- 
lence bonds). It is therefore appropriate to consider B{vi) 



as the "quantum dimer" operator to be used in place of 
the dimer density n x (vi) £ {0, 1} in the CDM. Because of 
the non-orthogonality of the valence-bond basis, D xx (r) 
is not, however, identical to the classical dimer-dimer cor- 
relation function. The two systems and their dimer cor- 
relation functions become identical in SU(iV) symmetric 
generalizations of the RVB when N — >• ooJ^ 

We will here show that D xx (r) for the standard S = 
1/2 SU(2) spins decays much slower than the classical 
correlator, as l/r a with a sa 1.20. These correlations, 
which are peaked at momenta q = (tt, 0) and q = (0, tt), 
correspond to critical fluctuations of a columnar valence- 
bond-solid (VBS). The exponent a < 2 in the RVB spin 
liquid corresponds to power-law divergent Bragg peaks, 
while in the CDM these peaks are only logarithmically 
divergent. As a consequence of the non-orthogonality of 
the valence-bond basis, the RVB is, thus, significantly 
closer to an ordered VBS state than is the CDM (or 
QDM). This result was first reported by us in a con- 
ference abstract^ and in an unpublished earlier version 
of this paper—, and was also found in independent par- 
allel work by Albuquerque and Aleti^ Here we provide 
further details on the dimer correlations and their signif- 
icance. 

We also study systems doped with two monomers and 
compute the distribution function of the monomer sep- 
aration. A well known result for the CDM is that the 
monomers are deconfined, with the distribution function 
M(r) decaying with the separation r as A{L)/r^ , where 
P = 1/2 and the prefactor A{L) decays with the sys- 
tem size L in such a way that the distribution is normal- 
ized for all L. For the RVB state, we find a more rapid 
power-law decay, with /3 ~ 0.83, which still corresponds 
to deconfined monomers. 

It is known that the dimer correlations of the CDM 
decay as 1/r 2 also in the presence of longer bipartite 
bonds (while non-bipartite bonds leads to a non-critical 
phase, with exponentially decaying correlations). As we 
will explain further below and in Appendix B, the ex- 
ponent a in this case does not correspond to these lead- 
ing correlations, however, but a subleading contribution 
decaying as l/r a with a > 2. This exponent and the 
monomer exponent j3 are non-universal, depending on 
details of the model (the fugacities corresponding to the 
longer bonds) ^ We also study here the RVB including 
longer bonds (the second bipartite bond, which connects 
fourth-nearest neighbors as considered previously in the 
CDM 1 -) and find that also in this case a and (3 change 
with the concentration of longer bonds. In contrast to 
the CDM, the leading dimer correlations arc always (at 
least for the range of parameters studied here) controlled 
by a, however, since a < 2 for the RVB. 



B. Height representation and topological sectors 

A key notion for relating the various results on the 
CDM and (we believe) the RVB model also, is that of 
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W=(0,0) W=(0,1) W=(0,3) 

FIG. 1: (Color online) Configurations in different winding 
number sectors, W = (W x ,W y ). Here W y is given by the 
number of bonds crossing the line drawn in the y-direction 
(since those bonds are at even y — shifting the bond configu- 
ration by one step in the ^-direction leads to W y —¥ —W y ). 
The last case is the unique configuration in its winding num- 
ber sector and constitutes the staggered state of the QDM.— 

"height model," or cquivalently a Z7(l) classical field the- 
ory. This means that all the long-wavelength behaviors of 
the system are captured by a coarse-grained scalar field 
h(r). The dimcr density operators and monomer defects 
can all be expressed in terms of h(r), and the weighting 
of its configurations is proportional to exp(— -Ftot); where 

Ftot = J d 2 r^K\\7h(r)\ 2 . (1.3) 

The height mapping for square-lattice dimers was in- 
troduced over twenty years agoj^lr— The use of such a 
mapping to explain correlation functions originated ear- 
lier (effectively for dimers on a honeycomb lattice) with 
Blotc, Hilhorst, and Nienhuis^ 

The key parameter in Eq. (jl.3|) is the dimensionlcss 
stiffness constant K. It can be shown that the expo- 
nents a and /3 measured in our simulations, as well as the 
coefficients of a "pinch-point" singularity in the dimcr- 
density structure factor, and also the ratios of the proba- 
bilities of different topological (winding number) sectors, 
are all functions purely of K. The details of the height- 
model construction underlying this result are given in 
Appendix [Bj It will be shown in Sec. [V] that all our 
measurements based on Monte Carlo simulations of the 
CDM and RVB consistently give the same value of K for 
a given model, demonstrating the validity of the height 
model. That is expected for the CDM, for which the 
height approach is well known; here we show that it is 
pertinent to the RVB as well. 

A related aspect of RVB states and the CDM is that 
their bond configurations on periodic lattices can be clas- 
sified according to a topological winding number^ We 
here define the winding number W = (W x ,W y ) as used 
in Ref. [2ft Drawing a path in the y direction, W y is the 
number of x-dimcrs crossed at even y minus the num- 
ber of such dimer crossing at odd y (see Fig. [IJ. An 
equivalent definition 6 - uses one of the W = single- 
domain states, such as the one in Fig.^a), as a reference 
state. As shown in Fig. [2jc), a direction can be assigned 
to loops of the transition graph so that each carries a 
"lattice flux"; if we call the net fluxes (§ x ,$ y ), then 
(W x ,W y ) ~ ($j,,$ x ) [or, depending on exactly which 




(a) (b) (c) 

FIG. 2: (Color online) (a) Reference state used here for defin- 
ing the winding number. The direction of the dimers is from 
sublattice B (open circle) to sublattice A (solid circle), (b) 
An arbitrary valence bond state, with dimers drawn in the 
opposite direction, from sublattice A to sublattice B. (c) The 
transition graph formed by the reference states in (a) and the 
arbitrary state in (b). The winding numbers correspond to 
the net fluxes (in units of the system length L) defined by 
traversing the loops formed along the arrows; here <f> x = 1 
and <& y = 0, or "!> = (1,0), which corresponds to winding 
number W — (0, 1) in the definition of Fig. \T\ 



reference state is used and how the y coordinates are as- 
signed, we could have (W x , W y ) = ($ y ,— & x ) — the signs 
are normally not important]. This definition can be di- 
rectly extended to systems with long dimers, by associat- 
ing that flux (which can have both x and y components, 
for cases where there are bonds not along the x or y 
axis) with a line connecting their endpoints. A third def- 
inition of the same winding number is (proportional to) 
the net height difference added up along a path crossing 
the system in the x or y direction, using the rules de- 
tailed in Appendix [Bj The possible winding values for an 
Lx L lattice are W x ,W y e {-L/2, -L/2 + 1, . . . , L/2}. 
The equal-weighted (CDM) ensemble is dominated by 
the winding number sector W = (W x ,W y ) = (0,0) [as 
follows from Vh = being the minimum of Eq. (|1.3I) ]. 

Recently, extended QDMs have been considered, with 
interaction terms that can drive the system into ground 
states with non-zero Vh in a sequence of commensurate 
locking transitions ! 25 ! 26 Quantum phase transitions in- 
volving these states are unusual, exhibiting aspects of 
deconfinement on a fractal curve of critical points (form- 
ing a Cantor set, which prompted the term "Cantor de- 
confinement" for this class of unconventional transitions). 
This motivates us to also study the CDM and RVB states 
in different winding number sectors, which (it turns out) 
also happens to be an effective probe of the states' topo- 
logical natures. In the case of the RVB, states defined 
within sectors of different winding numbers are not or- 
thogonal, but become orthogonal in the limit of the in- 
finite lattice (which we will here demonstrate explicitly 
based on simulations). 



C. Outline of the paper 

The outline of the rest of the paper is as follows: In 
Sec. |ll] we review the essential features of the valence 
bond basis that we use for the RVB-state calculations, 
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in particular how to extract spin correlations. The four- 
spin correlations are re-derived in detail in Appendix 
in an alternative way to a previous treatment of more 
general multi-spin interactions^! In Sec. IIIII we discuss 
Monte Carlo two-bond reconfiguration^ and loop-cluster 
algorithms for sampling the CDM and RVB states. We 
also discuss the winding numbers and issues related to 
sampling them cither grand-canonically (where there are 
some ergodicity issues in the case of the RVB) or canon- 
ically. In Sec. [TV] we present results for the standard 
case of only length- 1 dimers and valence bonds, as well 
as extended models with bonds of length \/5. In Sec. |V] 
the results are interpreted in terms of a height model. 
Detailed derivations of height model predictions are left 
to Appendix |B] In Sec. IVII we further characterize the 
nature of the critical VBS fluctuations in terms of the 
joint probability distribution of the order parameters for 
horizontal and vertical bond ordering. We conclude in 
Sec. IVIII with a brief summary and discussion. 

II. THE VALENCE BOND BASIS 

We work in the standard bipartite valence bond ba- 
sis, where a state of N (an even number of) spins on a 
bipartite lattice, 

N/2 

\ V ") = ^Nj4 lid - I I 2 ' 1 ) 

i=l 

is a product of singlets, where the first spin i of each 
singlet is on sublattice A and the second spin a(i) is on 
sublattice B. With the B sites also labeled as 1, ... , N/2, 
the set a(l), . . . , a(N/2) is a permutation of these num- 
bers and the label a = 1, . . . (A/2)! in \V a ) simply refers 
to all these permutations. The signs of the expansion 
coefficients of this state in the standard t, I s ~P m basis 
correspond to Marshall's sign rule for the ground state 
) °f a bipartite system^ i.e., 

Bign[* (^,...,^)] = (-ir^, (2.2) 

where is the number of J, spins on sublattice A. 

An amplitude-product state is a superposition of va- 
lence bond states, 

|*> = J2i>M, (2.3) 

a 

where the expansion coefficients are products of ampli- 
tudes h(r a> i) corresponding to the "shape" of the bonds 
(the bond lengths in the x and y direction in the case of 
a 2D system); 

N/2 

i>« = J! h{r ati ). (2.4) 
i=l 

Our main focus here will be on the extreme RVB state 
made up of only bonds of length 1 (one lattice constant), 




FIG. 3: (Color online) Two valence-bond states (left and 
right) in two dimensions and their transition graph formed 
by superimposing the two bond configurations (center). One 
of the spin configurations compatible with the transition 
graph is also shown, with open and solid circles for f and 
J, spins. Each loop has two such allowed staggered spin con- 
figuration, and the overlap of two valence-bond states is thus 
(Vp\V a ) = 2"°' 3 " JV/2 , here with the number of loops n Q(3 = 4 
and the number of spins N = 16. 

in which case the expansion coefficients ip a are the same 
for all configurations. We will also later study states 
including the bipartite bonds of length y/5 lattice con- 
stants, examples of which are seen in Fig. [3] The discus- 
sion here and in Sec. IIIII will be framed around generic 
bipartite amplitude-product states, with no restriction 
on the bond lengths. 

A. Transition graphs 

An important concept in the valence bond basis is the 
transition graph formed when the bond configurations of 
the two states are superimposed^^ This is illustrated in 
Fig- El The overlap (VglV^) between two valence-bond 
basis states can be simply expressed in terms of the num- 
ber n a p of loops in the transition graph. 

The easiest way to calculate the overlap is to go back 
to the standard basis of f and \, spins, so that 

(Vp\V a ) = ^££(-l)»^i+«^i x (2.5) 

si 

(Spi, . . . , Sp N \S al , . . . , S^ N ), 

where S a and S% denote spin configurations compatible 
with the bond configurations V a and Vp, i.e., those that 
have spins or 4-t on each bond. Terms with any oc- 
currence of S ai Sp t of course vanish, and the double 
sum, thus, simply counts the number of spin configu- 
rations common to the two bond configurations. Since 
the spins on each bond are antiparallel, the spins along a 
loop of alternating V a and Vp bonds (i.e., the loops in the 
transition graph) must alternate in a staggered, tltl ■ ■ •> 
pattern. There arc two such configurations for each loop. 
The total number of contributing spin configurations is 
therefore 2 n =< 3 , giving the overlap 

(Vp\V a ) = 2^~ N / 2 \ (2.6) 

which replaces the orthonormality condition (/3\a) = S a p 
for an orthonormal basis. For bond tilings V a = Vp, we 
have n a p — N/2 and the overlap equals unity. 
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In calculations with superpositions \ip) of valence-bond 
states, such as amplitude-product states, it is often not 
practical to normalize the states. It is convenient to write 
operator expectation values in the form 



where <f>ij is the staggered phase factor; 



(*|0|*> 



E a p^ a (Vp\0\V Q ) 



(2.7) 



Defining the weight W a p for the combined bond config- 
uration V a ,Vp and the normalized matrix element O a p 
according to 



W a p 

O a fS 



^ a {Vp\V a ) 

(Vp\Q\Va) 
(Vp\V a ) ' 



(2.8) 
(2.9) 



the expectation value takes the form appropriate for use 
with the Monte Carlo sampling methods that we will 
discuss below in Sec. Mil 



mom 



(2.10) 



The weight W a p, which is used in sampling the states 
in Monte Carlo simulations, is positive-definite when we 
consider wave functions satisfying Marshall's sign rule, 
i.e., the amplitudes h(r a> i) > in Eq. (|2.4[) . 

Like the overlap of the valence-bond states, the matrix 
elements of operators of interest can typically also be 
expressed in terms of the loops of the transition graph 
of the bond configuration V a ,Vp. We discuss spin and 
dimer correlations next. 



B. Correlation functions 

The standard spin-spin correlation function is most 
easily obtained by reintroducing the spins in the tran- 
sition graph, as illustrated in Fig. [3J We can then use 
the fact that 



s 3 \v a 



z(y p \szs z Av a 



(2.11) 



where the latter is diagonal and easy to compute in the z- 
spin basis. When summing over the allowed spin states, 
i.e., the two "orientations" of each loop (for a total of 
2™a/3 S pj n states), it is clear that S^S* averages to zero 
if i and j are in different loops, whereas for i,j in the 
same loop we get ±j(Vfs\V a ) , with the sign depending on 
whether the spins are in the same (+ sign) or different (— 
sign) sublattices. Introducing the notion for two 

spins in the same loop and for spins in different 

loops, we can write the matrix element ratio in Eq. (|2.10|) 
corresponding to the spin correlation function as 



(Vp\V a ) 



4 rij i 



(i,j)L, 



(2.12) 



-1, 
+1, 



for i,j on different sublattices, 
for i,j on the same sublattice. 



(2.13) 



While the loop-expression Eq. (|2.12[) for the simple spin- 
spin correlation function is well known,— ^ the general 
form of a four-spin correlation (of which the dimer-dimer 
correlator of interest here is a special case) was only de- 
rived recently^ In Appendix [X] we discuss this derivation 
in a slightly different way, which is less convenient when 
generalizing to higher-order correlators (which was also 
done in Ref. [27h , but more transparent in the case of the 
four-spin correlator. The resulting general formula for 
any non-zero four-spin matrix element is 



(^|(s fc -s ; )(s t -s J -)|y a ) 

(Vp\V a ) 

" (&-!$)^K> 

jg<i>ij<t>kl, 



(i,j) L (k,l) L , 
{i,k) L (j,l) L , 

(»)0i(j) k )L- 



(2.14) 



Here we have generalized the notation of Eq. (|2.12[) for 
how the sites are distributed among loops in a straight- 
forward way, with indices within the same parentheses 
belonging to the same loop. In the case of the single-loop 
contribution, (i,j,k,l)L, the term 5^ € {0,1} depends 
on the order of the four indices within the single loop, as 
specified in Eq. (|A9j) of Appendix [3] 



III. MONTE CARLO ALGORITHMS 

A simple but powerful Monte Carlo sampling algo- 
rithm for amplitude-product states based on reconfigu- 
ration of bond pairs was presented some times ago by 
Liang et al.^ who used this method to study the spin- 
spin correlations in amplitude-product states with sev- 
eral different forms of the amplitudes (exponentially or 
power-law decaying with the length of the bond) . A more 
efficient algorithm using loop updates was developed re- 
cently which operates in a combined basis of both valence 
bonds and spins.— The two-bond update, as well, can be 
made more efficient by working in this combined basis. 
Here we briefly review these two algorithms, and also dis- 
cuss the topological winding numbers that can be used 
to classify the bond configurations. 



A. Combined bond-spin basis 

Monte Carlo sampling of valence bonds involves mak- 
ing some change in the bra and ket bond configurations 
V a and Vp, and accepting or rejecting the update based 
on the change in the sampling weight Eq. (|2.8[) . accord- 
ing to some scheme satisfying detailed balance. Working 
with the standard non-orthogonal valence bond basis and 



using the Metropolis algorithm, we need to compute the 
weight ratio appearing in the acceptance probability 



P: 



accept 



• i 



(3.1) 



where the primes indicate the new states after some 
changes have been made in either bond configuration V a 
or V/3 (or both, but typically one would change only one 
state at a time). 

The weight ratio using Eq. ()2.8j) is 



0' 



W af 3 



lpalp/3 



(3.2) 



For an amplitude-product state, the ratio of the wave 
function coefficients is trivial, but computing the change 
n a i pi — n a p in the number of loops in the transition graph 
can be time consuming, as it involves tracing loops that 
can be long. 

The loops are typically long, O(N), if there is anti- 
ferromagnetic long-range orderJi That is not the case 
for the short-bond RVB states studied in this paper, but 
nevertheless it is more efficient to avoid the loop-counting 
step. That can simply be done by expressing each sin- 
glet in the standard basis of f and \. spins, and sampling 
these spin configurations in addition to the bond con- 
figurations [and since the spin basis is orthonormal, the 
sampled (non-zero weight) spin configurations must be 
the same in the bra and the ket]. That is, the configu- 
rations being sampled consist of a direct product of two 
valence bond patterns V a and Va, as well as one spin con- 
figuration Z a p compatible with both a and /3 (i.e. one ~f 
and one J. spin on each bond). Each loop in the transition 
graph must consist of an alternating string fit ••■ I anc h 
for every loop, there are two choices for this string. Thus, 
the ratio of the number of spin configurations is equal to 
the factor 2< n «'fi'-«» l !) i n Eq. $£2$. The Monte Carlo 
sampling of the spin configurations compatible with the 
bond configurations therefore automatically takes care of 
the factor 2™°< 3 in Eq. p. 21) . with no need to generate a 
transition graph or count loops. For more details of the 
arguments leading to this conclusion, see Ref. [l7l 



B. Monte Carlo sampling 




(a) 



(b) 



(c) 



FIG. 4: (Color online) Monte Carlo updates for the RVB 
state in the combined spin-bond basis. Open and solid circles 
represent f and 4- spins. In the basic moves (a) and (b), only 
one of the two two valence bond configurations is affected at 
a time, (a) A simple two-bond update. Choosing two sites 
on the same sublattice, the two bonds connected to them can 
be reconfigured in a unique way. If the spins are compatible 
with the f, J, singlet restriction, this update can be accepted, 
(b) Loop-cluster update. Choosing an arbitrary starting site 
(in this example in the left-upper corner) two defects (a site 
with no dimer or two dimers connected to it, both indicated 
with an x ) are generated by moving the end of the dimer on 
the initial site to another site which satisfies the bond-length 
constraint (here, in the extreme short-bond RVB, the length 
is always one) and the spin-singlet compatibility (anti-parallel 
spins on the bond) . The dimer that was previously connected 
to this site is then moved away from the double-bond defect 
to another site. This process continues until a bond returns 
back to "annihilate" the original empty-site defect, which here 
happens already after two bond moves [the last step in (b)]. In 
both (a) and (b) , we only show the bonds of the configuration 
involved in this update, (c) Monomer update. Monomers are 
shown as larger circles and must appear in the same locations 
in the state \V a ) and |Va), the bonds of both of which are 
shown here (as solid and dashed lines). In addition to the 
two-bond or loop update of the bonds, monomers can move 
to a site on the same sublattice by also moving a bond which 
is common to the two valence bond states. 



configurations. For either algorithm, updates are alter- 
nated between the ket and bra configurations, and there 
is an additional step for updating the spin configuration, 
where all the spins belonging to randomly chosen indi- 
vidual loops in the transition graph arc flipped. 



Here we outline the two different bond sampling al- 
gorithms that we used, each of which comes in a sim- 
ple version for the CDM, as well as a generalization for 
the combined spin-bond basis for the RVB amplitude- 
product states. In the case of the RVB, the spin config- 
urations also have to be updated. We also introduce a 
simple extension to sample states with monomers (empty 
sites). 

The two updating algorithms are summarized using 
simple examples with short bonds in Figs.2Ja,b), with (c) 
showing the extension needed for also sampling monomer 



1. Two-bond update 

For the two-bond update, as in Ref. [f| we choose two 
sites on the same sublattice (normally a next-nearest- 
neighbor site pair) and exchange their dimers in the 
unique way maintaining the A — B sublattice connectiv- 
ity, as shown in Fig. IDJa). The update can be accepted 
only if the spin configuration is compatible with the new 
bond structure, i.e., only antiparallel spins are connected 
by the bonds. In the case of the extreme short-bond 
RVB, an allowed new configuration is always accepted, 
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as the wave function ratio in Eq. (|3.2[) trivially equals 
one, whereas in general, when longer bonds are present, 
a ratio involving the amplitudes of two bonds has to be 
computed to determine the Metropolis acceptance rate 
Eq. (l3~Ti 

The algorithm for the CDM is simpler, as there is no 
spin state in that case. In the case of short bonds, an 
update of two bonds [flipping a pair of parallel bonds as 
in Fig.|4l[a)] is then always accepted, whereas in the pres- 
ence of longer bonds the acceptance probability involves 
the ratio of bond fugacities. We here consider only two 
bond lengths (nearest neighbor and fourth-nearest neigh- 
bor bonds, as shown in Fig. [3]), with fugacities Z\(i) = 1 
and Zi(i\ respectively, for bonds connected to site i 
(taken to be the sublattice A site, for definiteness). The 
partition function is then given by 



-CVM 



n 2 (C) 
2 



(3.3) 



c 



where 77.2(C) is the number of long bonds in configuration 
C. The acceptance probability for an update of bonds on 
sites i and j is 



Pxc.c = min 



Z 



'new (O^ncw (j ) 

Z id{i)Z id(j) 



(3.4) 



where "old" and "new" correspond to the length-index 1 
or 2 before and after the bond reconfiguration. 

For both the RVB and CDM, this algorithm keeps the 
system in a sector of fixed winding number, which we 
can take advantage of if we want to study properties in 
the individual sectors. Suitable starting configurations 
for different winding number sectors are shown in Fig. m 



2. Loop update 



that bonds have been moved on a closed loop of sites. A 
sweep of bond updates is defined as the construction of a 
fixed number of loops (determined during the equilibra- 
tion part of the simulation) which on average result in 
« TV moved bonds in both the ket and the bra state. 



3. Spin update 

After updating the bond configurations with one of 
the above algorithms, we update the spin configuration 
by flipping the spins of randomly selected loops of the 
transition graph (such as those in the middle graph of 
Figs. [3]), with probability 1/2 for each loop. All the loops 
have to be traversed, by moving between spins according 
to the bonds (which are stored in the computer as bidi- 
rectional links) , alternating between bonds in the bra and 
ket state. Each site visited is flagged and no new loops 
are started from already visited sites. The computational 
cost of a full sweep of such updates (visiting each site 
once) is 0(N). 



4- Monte Carlo sweep 

A sequence of bond updates in which 0{N) bonds are 
affected followed by a complete spin update constitutes 
one Monte Carlo sweep, which has a total computational 
cost O(N). Note that the sampling algorithm without 
the spins potentially costs up to TV 2 steps per sweep, 
since each two-bond update requires loop-traversals to 
check whether two loops are joined or a single loop is 
split r- and the loop length can then be up to O(N) (in a 
Neel state). The same issue pertains to loop updates in 
the pure valence-bond basis as well. 



If we want the system to wander among the different 
topological sectors, we instead use the loop-cluster up- 
date, which is a simple extension of a loop update for 
the CDM.— ^ It is also in general more efficient (ex- 
hibits shorter autocorrelation times) than the two-bond 
update for large size system. To start the loop update, 
we pick a site at random; in the example in Fig. |4|b) the 
top left site. We move the dimer connected to it, thus 
creating two defects in the system. We keep the starting 
site as a vacancy and move the original dimer of the now 
doubly occupied site to a new site, with certain probabil- 
ities satisfying detailed balance, and constrained by the 
spin configuration so that spins are opposite on every 
dimer. In the case of short bonds only, the probabili- 
ties are equal for the three new neighbor sites. For the 
general case where longer bonds are included, we refer 
to Ref. [TH for efficient choices of the probabilities. This 
update moves the doubly-occupied defect to a new site, 
which in Fig- HKb) is the lower-right site. We keep moving 
this defect using the above procedures, until it happens 
that the two defects annihilate each other, which means 



5. Sampling with monomers 

We will also be interested in the distribution of two 
monomers in the RVB states. In the case of the CDM, 
the distribution function of the monomer separation can 
be measured just by keeping track of the two defectS ) 12 i 30 
but in the RVB we have to explicitly introduce two 
monomers by removing both spins on a randomly chosen 
valence bond which is common to both the ket and bra 
bond configurations. Note that valence bond states with 
monomers are orthogonal unless the monomers are at the 
same locations in both states. We use the loop algorithm 
to sample the bond configuration space, and periodically 
we also move the monomers. Such a move can be done 
in combination with the move of a valence bond that is 
common to the two states, as shown in Fig. 0|c). This 
can always be accepted if there is no change in the bond 
length (one could also consider updates where a monomer 
moves and a bond length changes, which we do not do 
here). We update the position of two monomers in turn 
after each sweep of bond updates, when possible, and 
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measure the distribution probabilities M(r) as a func- 
tion of distance r between the two monomers. 

Note that if we assign spins to the monomer the situa- 
tion is different, due to the ovcrcompletcness of the basis. 
In a system with, e.g., two unpaired f spins, these two 
spins do not have to be located at the same sites in the ket 
and bra state — for a non-zero overlap it is only required 
that they are pairwise connected by valence bonds in the 
transition graph (which now contains two broken loops 
with open ends terminated by the unpaired spins). Such 
states with unpaired spins should be related to spinons, 2 
but we will not pursue studies of them here. Valence 
bond states including unpaired spins have recently been 
studied in different systemsi 31 i 32 



C. Winding numbers 

A two-bond update cannot bring the system from one 
topological winding number sector to another, while the 
loop update can. In the case of the RVB, there are 
winding numbers both for the bra and the ket state, 
and because of the non-orthogonality of the basis these 
winding numbers can be different. We denote the full 
winding number of a configuration in this case as W = 
(W" , Wy ; Wg , Wy 3 ) . In a grand canonical ensemble of all 
winding numbers, the sectors have different weight, which 
can be computed using Monte Carlo sampling with the 
loop updates simply by keeping track of the number of 
configurations generated in each sector. Results for such 
weights are presented below in Sec. IIV Al 

The loop algorithm for the CDM remains ergodic in 
the grand-canonical winding-number space even for very 
large systems, i.e., the loops can easily become very long 
and span the system. These long loops are related to de- 
confined monomers i^ 3 - The RVB simulations, in the case 
of short-bond states, in practice become stuck in some 
fixed winding-number sector for large L. However, the 
shortness of the RVB loops does not imply monomer 
confinement, as these loops are not directly related to 
states with monomers^ 3 - The loops for short-bond two- 
dimensional RVB states are typically very short (rarely 
exceeding 12 bonds in the case of the lcngth-1 bonds 
only) . This results in rather large error bars for computed 
quantities for L > 50, seen in grand-canonical results to 
be discussed further below. In practice, for large systems 
we will therefore study canonical ensembles in different 
fixed winding number sectors. Starting with a configu- 
ration initially prepared with a desired winding number 
(such as those illustrated in Fig. [I}, two-bond updates ex- 
plicitly conserve the winding number while loop updates 
in practice do as well, for large systems within reasonable 
simulation times. 



IV. RESULTS 

The ground state of the QDM at the RK point is the 
equal amplitude superposition of classical dimer states. 
The CDM can therefore give some insights into prop- 
erties of the RVB system as well, as long as the non- 
orthogonality of the valence-bond basis (i.e., the internal 
singlet structure of the valence bonds of the RVB) does 
not play an important role^ The quantitative validity 
of this approach is tested here by comparing the prop- 
erties of the CDM and the short-bond RVB state. We 
present the winding number distributions of both models 
m Sec. lRTAl then briefly discuss the standard spin corre- 
lation function of the RVB in Sec. lTVBl In Scc. lTTCl we 
study the four-spin VBS correlation function Eq. (jl.ll) of 
the RVB (which we also refer to as a dimer-dimer corre- 
lation function) and compare with analogous results for 
the well known dimer-dimer correlations of the CDM. 
In this section we consider the winding number sector 
W = (0, 0) and later, in Sec. IIVD1 discuss also cor- 
relations in systems with nonzero winding number. In 
Sec. HVEl we study the monomer distribution functions 
and in Sec. IIV Fl systems including the longer bonds. 

A. Sector probabilities 

We simulated the grand-canonical ensemble of wind- 
ing numbers, as explained in Sec. lIII Cl and accumulated 
the probabilities of several different sectors as shown in 
Fig. El for both the RVB and CDM, and for various 
system sizes L. The W = [(0, 0) for the CDM and 
(0, 0; 0, 0) for the RVB) sector is dominant in both cases, 
with the probabilities in the higher-W sectors decreas- 
ing rapidly. The probabilities of these \ow-W sectors 
clearly converge to L-indepcndent non-zero constants, 
rapidly with L for the CDM, and also for the diagonal 
(W a = W p ) sectors of the RVB (although the RVB data 
are much noisier for the large systems). By contrast, 
the probabilities of the off-diagonal sectors of the RVB, 
here exemplified by W = (0,1;0,0), decay exponentially 
to zero, which reflects the expectation that the states in 
different winding number sectors should become orthog- 
onal in the thermodynamic limits In the following, when 
considering winding number sectors of the RVB we will 
focus on the diagonal sectors and for simplicity denote 
the total winding number by W = (W x , W y ) in the same 
way as for the CDM. 

B. Spin correlations in the RVB state 

The spin-spin correlation function of the RVB has been 
studied before and is known to decay exponentially for 
a 2D system with short bonds (while a system with suf- 
ficiently slow decay of the probability of long bonds has 
long-range antiferromagnetic order) Here, we only 
comment briefly on the role of the winding number. For 
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FIG. 5: (Color online) Winding number probabilities obtained 
in simulations with the loop algorithms for the RVB and 
CDM (with only the shortest bonds, as in Fig. [4j. Results 
for several low-VK sectors of the CDM (lower panel) and RVB 
(upper panel) are shown versus the lattice size on a lin-log 
scale. In the RVB, the probability of the off-diagonal sector 
W — (0, 1;0, 0) vanishes exponentially with L, reflecting the 
orthogonality (when L — > oo) of states in different winding 
number sectors. 



unequal x and y winding numbers, W x ^ W y , the CDM 
and RVB systems clearly do not have the 90° rotational 
symmetry of the square lattice. We will investigate the 
directional dependence of the four-spin dimcr-dimer cor- 
relations below. Here, in Fig. [SJ we show results for the 
spin-spin correlations in two different winding number 
sectors. The correlations are always exponentially de- 
caying with distance, with a faster decay in the same 
direction as the one in which a non-zero winding number 
is imposed. 



C. Dimer Correlations 

In the CDM, the dimer-dimer correlation function 
D xx (r) is defined in the standard way using the bond 
occupation number n x (i) = 0,1 on the link of the lat- 
tice between site i and its neighbor at distance (1,0); 
D xx {^ij) = (nirij). The four-spin correlation function 
Eq. (jl.lj) of the RVB instead involves the loop estimator 
Eq. (TjOIl) ■ This reduces to the CDM form for SXJ(N) 
spins when N — > oo and the basis becomes orthogonal 
[in the representation of SU(./V) in which the factor 1/2 
in the off-diagonal matrix element in Eqs. (|A3[) and (|A4j) 
is replaced by see, Ref. HI for computations with 




FIG. 6: (Color online) Spin correlations versus lattice distance 
r in the short-bond RVB in the sector of winding numbers 
W = (0, 0) (top panel) and W = (0, 3L/7) (bottom panel) 
computed using L x L lattices with L = 48. Results are shown 
for the separation (x, y) taken along the two axis, (r, 0), (0, r), 
as well as on the diagonal, (r/v^, r /\/~2). 



such basis states]. For TV = 2, considered here, significant 
differences between the RVB and CDM can be expected. 

Since we are using periodic boundary conditions, the 
maximal separation to be used in the correlation func- 
tion is (L/2, L/2) on a L x L lattice. We first investi- 
gate the dominant part of the correlation function, which 
in the CDM is a mixture of a staggered component, at 
q = (ir, 7r) in reciprocal space, and columnar correlations, 
at q = (ir,0) and at (0,ir)J£ The asymptotic decay of 
these correlations can be accessed through the difference 
between the real-space correlations at two distances, e.g., 



D *xx{ x iV) = D xx (x,y) - D xx (x - l,y). 



(4.1) 



This quantity at the longest distance r = (L/2, L/2) is 
graphed versus L in Fig. [7] for both the RVB and the 
CDM in several fixed winding number sectors. 

For the CDM, the decay with L is consistent with the 
known ~ 1/r 2 decay of the dominant correlations. Apart 
from an overall prcfactor that depends on the winding 
number, there are only minor differences between the dif- 
ferent winding sectors for small systems. The dependence 
of the results on the winding number is stronger for the 
RVB, but, as expected, also here the exponent a in the 
power-law form 1 jr a becomes independent of W for large 
L (as long as the relative winding number W/L — > when 
L —> oo). Unlike the CDM, in this case the prefactor of 
the power-law form also converges as L — > oo, i.e., the 
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FIG. 7: (Color online) Dimer-dimer correlation function dif- 
ference Eq. (|4.1|l at the maximal distance versus the lattice 
size. The upper panel shows results for the quantum RVB in 
different topological sectors as well as in the grand canonical 
ensemble (including all winding number sectors, in which case 
the fluctuations between sectors becomes very slow for large 
systems, as reflected in the large error bar for L = 48). All 
correlations converge to the same power-law decay as system 
size increases. The power, based on the W = (0, 0) data for 
large L, is a = 1.191(6). The lower panel shows results for 
the CDM, which are consistent with ~ 1/r 2 (shown with the 
solid line) for all winding number sectors. 



correction to the prefactor decays as some power higher 
than a. 

In Fig. [7j wc also show results in the grand- 
canonical winding number ensemble, which, as discussed 
in Sec. lIII Cl suffers from problems with non-ergodic sam- 
pling for L > 50 (reflected in the large error bar for 
L = 48). For extracting the asymptotic form of the corre- 
lations, the W = (0, 0) sector is the best choice and gives 
D(r) oc l/r a with a = 1.191(6) for large systems. While 
the behavior is, thus, qualitatively similar to the CDM, 
the exponent differs considerably. The reduced value of 
the exponent can be interpreted as the RVB state being 
closer to an ordered VBS than might have been antici- 
pated based on the known CDM dimer correlations. 

There are two sources of differences between the cor- 
relations in the CDM and the RVB: the form of the esti- 
mator Eq. (|2.14p as well as the weighting of the bra and 
ket valence bond states with the loop factor 2™ Q < 3 for the 
RVB instead of the equal superposition of the individual 
bond configurations in the CDM. We have also measured 
the dimer correlations of the RVB in the same way as 
in the CDM, by just using the bond occupation num- 
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FIG. 8: (Color online) Fourier transform S(q) of the dimer- 
dimer correlation function D xx (r) for systems of size L — 
32. The squares represent the full reciprocal space q x ,q y £ 
[0, 27r]. Results in winding number sectors W — (0,0), W = 
(0, 1), and W = (0, 8) are shown for the RVB (left) and CDM 
(right). The location of the broad ("incommensurate") peak 
in both cases is Q = (n, 2nW y /L). The sharp peak at (n, n) 
is due to a nonzero average staggered dimer order induced by 
a nonzero winding number. This peak has been removed in 
the graphs W = (0, 8) in order to make the other features of 
the correlations better visible. The height of the peaks as a 
function of the system size is analyzed in Fig. [9] 



bers in the bra and the ket states (but with the correctly 
weighted sampling of the RVB). We find the same ex- 
ponent a m 1.20 as above, which shows that the source 
of the different power-law is only the different weight- 
ing of the states. This could also have been anticipated 
based on the fact that the spin-spin correlation function 
of the RVB is exponentially decaying, which translates 
into short loops in the transition graph^ The loop esti- 
mator Eq. (|2.14[) of the four-spin dimer correlation func- 
tion is therefore still local and cannot change a power 
law. 

The Fourier transform of the full dimer-dimer cor- 
relation function D xx (r) is the structure factor S(q). 
This quantity gives a more detailed picture of the long- 
distance behavior of the dominant correlations. Rcprc- 
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FIG. 9: (Color online) Peak values of the dimer structure fac- 
tor, where Q = (n, 2nW y /L), versus the system size in sectors 
with different winding W y . The modified definition S (Q) for 
the RVB is given in Eq. (|4,2|l . Note the different j/-axis scales 
used for the two models (logarithmic for the RVB and lin- 
ear for the CDM). In the CDM (lower panel) the behavior is 
consistent with a log divergence (as shown with fitted lines) 
for small winding numbers, but for larger W it appears that 
the behavior is instead governed by a power law (which then 
may be the case for all W x /L > for sufficiently large sys- 
tems). The curve through the W = (0,3L/7) data shows 
5(Q) oc L ' 48 . In the RVB (upper panel) the exponent of the 
power-law divergence decreases slightly with increasing wind- 
ing number. The legends with (S) correspond to the peak 
values of the full structure factor S(Q). 



scntativc results for the 5(q) for L = 32 systems in three 
different winding number sectors (0, W y ) are shown in 
Fig. [5] In this section we focus on the W = (0, 0) sec- 
tor and leave discussions of nonzero winding numbers to 
Sec. ITVDl The "bow-tie" feature seen for W = (0, 0) 
in the CDM is well known and understood based on 
the mapping of the system to a height model (see Ap- 
pendix [B]) . The system has two kinds of power- law cor- 
relations: an effectively dipolar kind, which is responsible 
for the "pinch-point" singularity at q — (tt, tt) (see Sec. 
IB 3[) . and a "critical" kind with variable exponents, which 
leads to a broad peak at Q = (tt, 0) diverging logarith- 
mically with the system size, as shown in the lower panel 
of Fig. [9] In the RVB the peak is much sharper and 
diverges faster, as a power law (as shown in the upper 
panel of Fig. [5]) on account of the real-space form 1 jr a 
with aw 1.2 < 2 of the dimer correlation function. 



When the Fourier transform S(q) is computed post- 
simulation based on all computed real-space correlations, 
the measurements in the simulations are expensive, re- 
quiring 0(N 2 ) operations to take full advantage of spa- 
tial averaging. In the CDM, we can instead easily just 
compute S'(Q) at the single wave- vector Q directly in 
the simulations at a much lower cost of O(N) to access 
larger system sizes. In the RVB, this speed-up is not pos- 
sible, however, because we are there really measuring a 
four-spin correlation function that cannot be simply ex- 
pressed as a product of two-spin correlators, as discussed 
in Appendix A, and there is no obvious way of avoiding 
the 0(N 2 ) scaling of this measurement. 

In order to have a similar quantity, which scales with 
the system size in the same way as S'(Q) but for which the 
measurements require only O(N) operations, we define a 
modified structure factor S"(Q) for the RVB as 

S'(Q) = (B* X (Q)B X (Q)) (4.2) 

where B X (Q) is the Fourier transform of the spin-spin 
correlator matrix element (Vg|(S,; • Sj)|V^) for an indi- 
vidual configuration in the RVB simulation (i.e., obtained 
from a transition graph, which gives values € {—3/4,0} 
for each nearest- neighbor bond on the lattice). This def- 
inition of the peak value differs from the full Fourier 
transform S*(Q) of the four-spin dimer correlator D(r), 
essentially because it does not contain any information 
on the order of the site indices in the matrix element 
(Vg|(Sfc • S;)(S ? ; • Sj)\V a ), which plays a role in the 
transition-graph two-loop estimator of the dimer correla- 
tion function (as discussed in Appendix A). In particular, 
the modified quantity misses certain negative contribu- 
tions arising in some cases where all four indices belong 
to the same loop [see Eq. (|2.14[) ]. Therefore, we expect 
S"(Q) > S(Q), which is also confirmed by results for 
both quantities in small systems, as shown in the upper 
panel of Fig. |H1 The form of the power-law divergence is 
the same, however. 

Overall, there is significant directional dependence in 
the dimer correlations, but for W = (0, 0) the RVB re- 
sults in Fig. [5] confirm that the peak at (tt, 0) (corre- 
sponding to columnar-modulated correlations) is suffi- 
ciently isotropic for the size dependence of the Fourier 
peak to be directly related to the exponent of the power- 
law decay 1 /r a found above for the real space correlation 
(and, it should be pointed out, the exponent a also comes 
out consistently to the same value when extracted in dif- 
ferent directions in real space). 

With S"(Q) diverging with the system size L as L ctQ , 
we expect cxq ~ 2 — a and the data confirm this. For 
instance, the W = (0,0) data in the upper panel of 
Fig. 02 was fitted to a function f(L) = b Q L a Q + b 2 L a2 , 
where a 2 < cxq (and typically also a 2 < 0) and this cor- 
rection term is added in order to include data for the 
full range of systsem sizes. By using this form we ob- 
tained olq = 0.800(2), which is in good agreement with 
a = 1.191(6) but with a smaller error bar. Our best es- 
timate for the exponent is, thus, a = 1.200(2). Here the 
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FIG. 10: (Color online) Dimer correlation differences versus 
x [where the separation r = (x, x)] along the diagonal lattice 
direction for systems of different size. The winding number 
is W = (0,2), and therefore two phase shifts are seen (corre- 
sponding to a total of four domains). Note that the overall 
magnitude of these correlations is much larger in the RVB 
(upper panel) than in the CDM (lower panel). 



L/2 



RVB W=(0,0) 



L/2 



L/2 



RVB W= (0,1) 



L/2 



error bar is purely statistical and there may still be some 
systematical errors present as well (likely of the same or- 
der), arising from neglected higher-order corrections. 



D. Correlations with nonzero winding number 

In a background of nonzero winding number, all dimer- 
dimer correlations should become modulated by the fac- 
tor cos(£Q -r), as derived using the height-model formal- 
ism in Appendix [B] and shown explicitly as Eq. (|B16I) , 
where SQ = 2ir(W x , W y )/L. Such a modulation is visi- 
ble in the real-space dimer correlation function, as shown 
in Fig. [TU] for D* x (r) along the diagonal lattice direc- 
tion, r = (x,x), for systems of different size with wind- 
ing number W = (0,2). This implies that when r is 
followed along the [1, ±1] direction through an entire pe- 
riod, 2(W X ± W v ) nodes of D xx (r) are crossed; indeed, 
Fig.[TU]foi' W = (0, 2) shows two changes of sign between 
x = and L/2, in both the CDM and the RVB cases. 

The correlation function D xx (x, y) in the full 2D space 
is shown for the RVB in Fig. [TTJ where an overall back- 
ground constant representing D(r — > oo) has been sub- 
tracted from D(r) and the remainder has been multiplied 
by r a to make the modulations visible. An over-all non- 
decaying staggered contribution present when W =/= 
has also been subtracted (see further discussion of this 
below and in Fig. [5]). The color coding shows positive 
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FIG. 11: (Color online) Correlation patterns obtained from 
the dimer correlator D xx (x,y) by subtracting a constant and 
dividing the result by the leading power-law form r~ a (with 
a = 1.2 for the RVB and a = 2 for the CDM). The (vr, tt) con- 
tribution was also removed for the W =^ sectors (by going 
to Fourier space as in Fig. [8j) . Black and red (gray) bars rep- 
resent positive and negative values (i.e., stronger and weaker 
dimer correlations), respectively. In the W — (0,0) sector, a 
dominant columnar pattern is visible, while in the W = (0, 1) 
sector the correlations shift from weak-strong weak-strong to 
strong-weak strong-weak over a window of distances oc L, cor- 
responding to two nodal lines as stated in text. The origin is 
at lower left corner, and one quadrant (L/2 x L/2) is shown 
of the possible separations. In the W — (1, 1) sector, corre- 
lations shift twice in a row, corresponding to the presence of 
two pairs of nodal lines. 
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and negative correlations, and the width of bars repre- 
sent the magnitude of the correlations. In the winding 
number W = (0, 0) sector, the positive and negative val- 
ues alternate in rows, showing that the overall dominant 
correlations are of columnar type. In the W = (0, 1) sec- 
tor, a phase shift occurring around at y = L/2 is clear. 
The region over which the shift takes place is itself of size 
O(L), as expected since the amplitude is modulated pro- 
portional to a sine wave (which can be considered as a 
highly fluctuating critical delocalized domain wall). The 
results for the W = (1,1) sector confirm the existence of 
two such delocalized nodes along the diagonal direction. 
A similar pattern of phase shifts in the correlation func- 
tion is seen in the CDM well, but is much weaker 
because of the significantly faster decaying correlations 
(as is also clear in Fig. [TU|) . 

To our knowledge, these correlations in sectors of fixed 
non-zero winding number have not been studied in de- 
tail previously (but were pointed out also in the parallel 
work by Albuquerque and A\e^2Q) . In Appendix [B] we 
extend the height-model approach to this case as well (in 
Sec. IB 7[) . Here we only briefly discuss some of the main 
features, with the aim of comparing the RVB and CDM 
systems. 

Turning back to the Fourier space plot, Fig. [5J it in- 
cludes representative results for the structure factor in 
three different winding number sectors (0, W y ). Once the 
winding number is non-zero, it is clear that there is, for 
both models, a 6- function peak in 5(q) at (ir, 7r), reflect- 
ing a non-zero static staggered order parameter. Since 
this peak grows in proportion to the winding number, we 
have subtracted it off in some cases in Fig. [3] to make the 
other features better visible. 

There are two notable features of these results, for both 
the RVB and CDM: (i) the pinch-point remains at (ir, ir) 
and (ii) the singularity at (it, 0) present for W y = is 
offset to Q = (tt, 2irW y /L), which when L — > oo can be 
considered as an incommensurate peak at Q = (tt, w), 
w G [0,7r]. This is exactly as expected from Eq. (|B17|) 
obtained within the height-model representation in Ap- 
pendix [Bj Figure O shows the system size dependence 
of the singular peak for different large winding numbers 
W y ex L. These features have been qualitatively ex- 
pected in the case of the CDM based on several previous 
work a 26 ' 36 ' 37 (as outlined in Appendix [Bj) . but they are 
still interesting to study quantitatively and to elucidate 
the similarities and differences between the CDM and 
RVB. It is already clear from Fig. [8] that the divergence 
of the incommensurate peaks is much stronger for the 
RVB than the CDM, which is anticipated based on our 
result for the slow real-space decay of the dimer-dimer 
correlations in the RVB. 

For non-zero winding number, the correlations become 
significantly anisotropic, but we have not attempted to 
study their full functional form in real space or Fourier 
space. The exponent governing the asymptotic power- 
law decay is, however, expected to be direction indepen- 
dent, as discussed in Appendix [Bj The results in Fig. [5] 
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FIG. 12: (Color online) Monomer distribution function in the 
RVB state on an L = 512 lattice. The straight line is a fit to 
the power-law form l/r fi with /3 = 0.830(9). 

indicate that S(Q) has the form L aQ , with a weak de- 
pendence of the exponent <xq on the location of the peak 
(i.e., the winding number), also as expected based on the 
height- model results in Appendix IB 81 

The incommensurate peak of the CDM was discussed 
by Fradkin et al.,— who pointed out a set of critical 
points in extended QDMs with more complicated diago- 
nal and off-diagonal terms than the standard RK nearest- 
neighbor bond-pair interactions. The critical points ex- 
tend from the conventional RK point at zero winding 
number, forming a complex fractal curve with devil's 
staircase features (forming a Cantor set). This critical 
curve separates a staggered dimer phase from one with 
a complex bond pattern with a large unit cell, which de- 
pends on the winding number. Similar transitions with 
a series of different VBS phases were studied in Ref. Hf| 
Our CDM results in Fig.[9]for large winding numbers sug- 
gest that the incommensurate peak may become power- 
law divergent (i.e., stronger than the logarithmic diver- 
gence obtaining at zero winding number). This is seen 
most clearly in the W = (0,3L/7) graph, where it is 
clear that the divergence with L is faster than logarith- 
mic. A power-law fit, L aQ with o.q = 0.48(3) describes 
the data well. This is expected in the height scenario, 
since a nonzero background W/L changes the effective 
stiffness to K' as given by Eq. (|B25[) . The exponent a 
of real-space correlations accordingly changes from 2 and 
consequently the integral of l/r a (the structure factor) 
should diverge faster than logarithmically. 

E. Monomer distribution 

Monomers are expected to be deconfincd in RVB 
statesji which provides an intuitive picture of spin-charge 
separation. Here we will study two monomers in the 
RVB. It should be noted, however, that these monomers 
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FIG. 13: (Color online) Dimer-dimer correlation function dif- 
ference D xx (r) for RVB systems in the W = (0, 0) sector with 
different fugacities Z% of long (fourth-neighbor) bonds (with 
the short-bond fugacity Z\ = 1). The decay exponents grows 
with the long-bond fugacity. The values are given in Tabled] 




r 



FIG. 14: (Color online) Monomer distribution function M(r) 
for RVB states with a small fraction of fourth-neighbor bonds 
on a lattice of size L = 256. The straight lines are fits giving 
deconfinement exponents which decrease with increasing long- 
bond fugacity. The exponents are listed in Table [I] 



are bosonic, and hence the results cannot be directly re- 
lated to a hole-doped RVB spin liquid. In that case the 
monomers should be fermions and, as discussed, e.g., in 
Ref. i, the sign rule we use here for the valence bonds 
would have to be replaced by more complex signs. It is 
nevertheless interesting to compare the monomer-doped 
RVB and CDM systems considered as different statistical 
mechanical systems. 

The monomer-monomer distribution function of the 
CDM is defined using the monomer density m(r,) = 0, 1; 



M{v l3 ) 



{m{Ti)m{Tj)) 
(m(r i )m(r l + x)) 



(4.3) 



where the normalization with the correlation at distance 
r = 1 is a convention which makes it easy to compare 
results for different system sizes (i.e., results for fixed 
r converge to a non-zero number with increasing size, 
even if the monomers are deconfincd) . It is knowni^ that 
this function for the short-bond CDM decays as M(r) oc 
1/r" with /3 = 1/2. This slow decay reflects monomer 
deconfinement, i.e., the function (m(rj)m(rj)) without 
the normalization in Eq. (|4.3|) decays to zero for fixed 
when K — > oc. We use exactly the same definition of 
M(r) for the RVB, applying the procedures discussed in 
Sec. IIIII to sample monomer configurations (while in the 
CDM the loop algorithm for the bond sampling without 
monomers gives the monomer distribution function as a 
by-produc t 12 i 30 ). Note that the winding number is not 
well defined in the presence of monomers, since they are 
associated with "broken loops" in the transition graph in 

Fig.n 

The exponent (3 = 1/2 for the CDM has been con- 
firmed previously in Monte Carlo simulations on large 
lattices^ Figure [T2l shows our results for the RVB, us- 
ing a system of size L = 512 (for which the results for 



moderate separation of the monomer are sufficiently con- 
verged to extract the decay exponent). We find that 
the exponent (3 sa 0.83 is significantly larger than in the 
CDM. The monomers are, thus, more strongly correlated 
to each other than in the CDM, but still deconfined. Note 
that in a long-range ordered VBS one would expect the 
monomers to be confined. 



F. Including longer bonds 

As the next step after investigating the extreme short- 
bond RVB, it is natural to think about the role of 
the longer bond in spin liquids and the classical dimcr 
model. In the case of the CDM, introducing bonds 
between next-nearest neighbors on the square lattice 
leads to exponentially decaying dimer correlations and 
monomer confinement,— as on a triangular lattice with 
only nearest-neighbor bonds However, with only bi- 
partite bonds, the behavior is qualitatively similar to 
the short-bond model (as long as the fugacity for longer 
bonds decays sufficiently rapidly with the length of the 
bonds) The dimer correlations decay as l/r a with 
a = 2 not changing as longer bonds are introduced, but 
the monomer exponent a decreases from 1/2. 

In the RVB, Marshall's sign rule cannot be applied 
if non-bipartite (frustrated) bonds are introduced. Due 
to the non-orthogonality of the basis, there is, regard- 
less of how signs beyond some simple Marshall rule are 
introduced, a sign problem in the Monte Carlo bonds 
sampling (due to non-positive definitcness of the state 
overlaps). We here study the effects of bipartite va- 
lence bonds connecting fourth- nearest neighbors, i.e., of 
"shape" (x, y) — (2, 1) and all symmetry-related shapes, 
as was done previously for the CDMJ^ We use small 
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TABLE I: Dimer-dimer and monomer exponents obtained for 
the CDM and RVB systems at different fugacities Z 2 for the 
next-shortest bonds (of length y/E). 



Model 


z 2 


a 


/3 


CDM 


o 


1.98(1) 


0.4996(5) 


CDM 


e" 4 


2.17(2) 


0.447(2) 


CDM 


e" 3 


2.44(8) 


0.392(1) 


CDM 


e' 2 


2.7(2) 


0.302(1) 


RVB 





1.191(6) 


0.830(9) 


RVB 


e" 4 


1.255(5) 


0.775(5) 


RVB 


e" 3 


1.377(10) 


0.707(5) 


RVB 


e" 2 


1.676(12) 


0.563(6) 



fugacities Z2 = e~ 2 , Z2 = e~ 3 and Z2 = e~ 4 for the 
longer dimers and Z± = 1 for the short bond connect- 
ing nearest neighbors. In the RVB, since we work with 
the amplitude product states [Eq. (|2.4p j. we just use the 
"fugacities" as another notation for the RVB amplitudes; 
h(r = 1) = Z 1 = 1, h(r = s/E) = Z 2 . 

Spin correlations have been previously studied in the 
presence of long bonds, including exponential and power- 
law decays of the length-dependent fugacities^^ Here 
we again focus on the dimer-dimer correlations and 
monomer distribution function. 

The exponent of the dimer-dimer correlations changes 
with the fugacity of long bonds, as shown in Fig, [inland 
Table |TJ The change can be seen even more obviously in 
higher winding number sectors (not shown in the figure). 
Note also that the spin correlations increase when longer 
bond are introduced£i2i Fig. PHI shows the monomer dis- 
tribution M(r) as defined in Eq. (|4.3j) . Similar to the 
CDMji^ the confinement exponent changes with fugacity 
of long bonds. The higher the fugacity of long bonds, the 
lower is the monomer deconfinement exponent. 

V. HEIGHT MODEL INTERPRETATION 

All of the numerical results found in these simulations 
can be compared with results obtained in the framework 
of the "height model" introduced in Sec. II Bl and elabo- 
rated in appendix [Bj According to that description, each 
of the following can be written as a function of a single 
parameter, the height stiffness K: 

(1) The sector probabilities P(W x ,W y ) presented in 
Fig. El 

(2) The exponent a of critical dimer correlations, in- 
ferred from the i-dcpcndcncc of the structure fac- 
tor at Q = (7r, 0) [the peak- value at winding num- 
ber W = (0,0) as shown in Fig. [9], and also from 
the L dependence of these same correlations at 
r = (L/2,L/2) in real space, as plotted in Fig. [7] 

(3) The decay exponent /3 of the monomer distribution 
function M(r) as presented in Fig. [T2l 



TABLE II: Stiffness parameter Kp in the infinite CDM and 
RVB systems inferred from the winding-number sector prob- 
abilities (from data in Fig. [5j according to Eq. (15.11) . 





CDM 




RVB 


(W X ,Wy) 


P(W X ,Wy) 


K P 


P{W X ,Wy) K P 


(0,0) 


0.49625(4) 




0.764(5) 


(1,0) 


0.10321(3) 


0.19628(3) 


0.057(2) 0.325(5) 


(1,1) 


0.02146(1) 


0.19629(3) 


0.0043(5) 0.324(7) 


(2,0) 


0.000925(2) 


0.19642(8) 





(4) The coefficient of the "pinch-point" singularity in 
the structure factor S(q) as shown in Fig. [8] 

We can use these relations to reduce the different results 
to independent estimates of the stiffness, which we call 
Kp, K a , Kp, and K$, from these respective measure- 
ments. The agreement (to be demonstrated below) of 
these is powerful evidence that a height-like field the- 
ory underlies the RVB state. That is well-known to be 
true for the CDM state, but the extension to the RVB is 
non-trivial, due to the configuration space here consisting 
of two bond configurations weighted by their transition- 
graph loops, as discussed in Sec. [TTJ Indeed, we have not 
derived the height-model representation explicitly for the 
RVB. We will make some comments on the feasibility of 
actually deriving the effective model below. 



A. Four ways to extract stiffness 

We now run through the ways in which we get four in- 
dependent measurements of the height stiffness K. CDM 
results are presented in parallel to the RVB results, firstly 
to check the systematic errors in our fitting procedures 
against exactly known results, and secondly to emphasize 
the similar behaviors. 



1. Sector probabilities 

Table [XT] gathers together the numerical sector prob- 
abilities from the data sets in Fig. [5j As seen in the 
figure, the smaller sizes show noticeable finite- L correc- 
tions, which are expected to be 0(1/ L 2 ) due to the quar- 
tic correction Eq. (|B19[) to the free energy density. The 
larger sizes show larger statistical errors particularly for 
the RVB explained in Sec. 1111 CI In order to 

partially account for finite-!/ corrections of leading order 
and higher, which we need to extract the probabilities at 
L — > 00 with relatively smaller statistical fluctuations by 
using a large set of lattice sizes, we use suitable polyno- 
mial fitting functions (some times without linear term) 
to extrapolate values in the thermodynamic limit. 

According to Eq. (|B15I) . we expect P(W x ,W y ) oc 
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0.230 



^ 0.225 



0.220 



0.215 



0.1167 



0.1164 



0.1161 



0.1158 




FIG. 15: (Color online) K P value calculated in the W = (0, 1) 
sector according to Eq. (|5.1|l for systems with fugacity Z2 = 
e~ 2 for long bonds and different lattice sizes. RVB and CDM 
results are shown in the upper and lower panel, respectively, 
as a function of the inverse system size 1/L. The curves are 
second-order polynomial fits, not including the linear term. 



cxp[-8K(W x + Wy)], and thus, we define 
ln[P(W x ,W y )/P(0,0)] 



Kp = - 



8{W% 



(5.1) 



This expression clearly gives consistent results for every 
pair (W x ,W y ), for either model as shown in Table HU 
The Kp values in this table are calculated directly from 
the corresponding sector probabilities presented next to 
them. The Kp values included in Table Mil are taken 
from the W = (0,1) sector, as that has the smallest er- 
ror bars (and also should be the best in terms of origi- 
nating from a weak "tilt" field). As indicated by Fig. H"5l 
the Kp value does not depend much on system size L 
for L larger than sa 50. Therefore, in order to obtain 
smaller statistical errors, we presented Kp in Table IIIII 
with the same method described above for extrapolating 
winding sector probabilities in the thermodynamic limit. 
As an example, polynomial fitting functions are shown in 

Fig.cna 



TABLE III: Stiffness estimates obtained from the four kinds 
of measurements discussed in the text; Z2 is the fugacity for 
dimers of length yE. 



Model 


Z2 


K P 


K a 


Kp 


Ks 


CDM 





0.19628(4) 


0.198(1) 


0.1962(2) 


0.1959(7) 


CDM 


e" 4 


0.17547(4) 


0.182(2) 


0.1755(8) 


0.1794(3) 


CDM 


e" 3 


0.15065(6) 


0.161(5) 


0.1539(4) 


0.1582(4) 


CDM 


e~ 2 


0.11638(3) 


0.14(1) 


0.1186(4) 


0.1234(1) 


RVB 





0.323(5) 


0.330(2) 


0.326(4) 


0.3242(4) 


RVB 


e" 4 


0.3067(8) 


0.313(1) 


0.304(2) 


0.3081(2) 


RVB 


e" 3 


0.2774(5) 


0.285(2) 


0.278(2) 


0.277(1) 


RVB 


e" 2 


0.2258(1) 


0.234(2) 


0.221(2) 


0.22619(2) 



The values of a summarized in Table U could in princi- 
ple all be obtained by fitting the size dependence of the 
peak- value S(Q) of the dimer structure factor, i.e., ac- 
cording to the peak-height analysis illustrated in Fig. [9] 
in the case of the RVBs. However, this approach requires 
a very significant computational effort for large lattices. 
We therefore use an easier but still reasonably accurate 
way to extract a, by fitting the real-space long-distance 
dimer correlator D* x (L/2, L/2) as in Fig.[T3]by a power- 
law [as expected according to Eq. (|B7[) ] . For non-zero 
Z2 cases in the CDM, this approach does not work well, 
however, because a increases with the fugacity, becoming 
larger than 2, and therefore the critical term is overshad- 
owed by the stronger dipolar term (which always decays 
as 1/r 2 ; sec Sec. IB 4[) and is hard to detect. In contrast, 
in the RVB a < 2 always and the critical term is domi- 
nant. A better way to find a in the CDM is to extract 
values by a fit of \D* x (x, x)\ (along the diagonal axis) for 
a range of distances x on a large lattice, since the dipolar 
term vanishes on this axis. The corresponding K a values 
are listed in Table IIIII 



3. Monomer pair distribution correlations 



We have [Eq. (|B~TT|) in Appendix El] that /3 = SK/tt; 
hence we define 



Kn = — — . 



(5.3) 



This quantity extracted from the exponents listed in 
Table U where the values originate from fits to 
the r-dependence of the monomer distribution function 
(Fig. [H in the case of the RVBs), is listed in Table UTTJ 



2. Critical dimer correlations 



4- Coefficient of the pinch-point in <S(q) 



We have [Eq. (|B8} in Appendix IFi] that a = tt/8K; 
hence we define 



K n 



- 

8a' 



(5.2) 



At Q = (jr, 7r), there is a pinch-point singularity of 
the dimer structure factor for cc-oriented dimers, 5(q), 
meaning that there is no divergence, but the limiting 
value at Q depends on the direction of the ray along 
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which it is approached. The coefficient of this ky/{k^. + 
ky) singularity is 1/K according to Eq. (|B5[) . so we can 
do a simple fit and call the result K$- Of course, the 
actual dependence on q = Q + k must have additions of 
higher order in k, since S(q L ) is periodic in the Brillouin 
zone. Therefore, only a small domain around Q should 
be used in the fit, but it may be advantageous to use 
more than the wave-vectors immediately adjacent to Q, 
as one can then extrapolate to Q and eliminate most of 
the unwanted additions. Of the four methods, this one 
is closest to direct measurement of the height Fourier 
spectrum (|/i(k)| 2 ), which was the best method to extract 
stiffness constants from simulations of height model a 37 ' 38 
or random-tiling quasicrystal o 39 ^ 40 . 

In the RVB case, some additional steps are necessary, 
because we do not construct a height function and do not 
really even have a dimer configuration (recall that the 
contributions to the wave function from different dimer 
configurations are non-orthogonal and the simulations 
sample pairs of dimer configurations). We only have the 
correlations D xx of an operator that has some projection 
onto a dimer-like variable as well as other contributions. 
This has two consequences for 5(q). The first is that the 
"other contributions" contribute a constant background 
on top of the pinch-point singularity, which does not van- 
ish even along the line k y = 0. That can in principle be 
remedied by fitting and subtracting off the constant addi- 
tion, but unless the lattice is very large such a procedure 
will not be perfect. In our fits carried out here, we simply 
use the value of the point that is next to the pinch point 
along line k y = as our constant addition. 

The second consequence of the lack of a formal height 
model is that the measured 5(q) is a multiple of the as- 
sumed dimer structure factor by an unknown coefficient 
c|. Fortunately, we can calibrate c| using the sectors 
with nonzero winding numbers, since the (5-function peak 
at Q in those cases (after subtracting the constant back- 
ground) is proportional to c| times (W 2 + W£) times 
known constants, allowing us to infer c| w 0.56. From 
this value we can extract a normalized S(q) and, finally, 
find the pinch-point coefficient we call 1/K$- This esti- 
mate of Ks was computed for several system sizes and 
then extrapolated to L = oo by fitting functions f(L) = 
a +a 1 /L 2 +a 2 /L 3 for the RVB and f(L) = a +ai/L 2 for 
the CDM (i.e., with both forms not including the linear 
term) . The results are given in Table Mil 

B. Summary of the stiffness estimates 

Table Hill collects all four estimates of K , with their 
statistical errors (one standard deviation). The fugacity 
Zi for long dimers specifies a family of RVB models and 
one of CDM models, with different exponents. Note that 
K according to our convention is 7r/8 times K as used 
previously in Ref. [H. 

The respective estimates for the stiffness constant for 
a given case typically agree to within a few error bars. In 



some cases the deviations are larger than expected purely 
based on statistics. This is not unexpected, since the cor- 
relation functions we have analyzed are also affected by 
corrections to the leading forms we have used. Note that 
Ks for the CDM with long dimers are systematically too 
large (the only really significant disagreement); and A'5 
for the RVB with long dimers appears to be slightly too 
large as well. Here the background contributions which 
may not be perfectly subtracted off in our procedure, 
may be to blame. 

The results for the CDM can be compared with the 
exact value A'cdm = n/16 w 0.19635, with which all 
K estimates in Table IIIII agree to within 2 error bars 
or less. As another test, we calculated Kp for the CDM 
with long bonds only (i.e., fugacitics Z2 = 1 and Z\ =0). 
The resulting value implies an exponent for the monomer 
correlations of j3 = 0.11092(6), which agrees (within 1.5 
error bars) with a previous obtained using a different 
analysis of the monomer distribution function (and where 
it was conjectures that fi = 1/9) «i£ 

The good agreement between four different stiffness es- 
timates provides strong evidence of an underlying height 
model description of the RVBs. The plausibility of the 
height-model approach for the RVB is partially motivated 
by the fact that the RVB and CDM coincide for SU(A) 
spins when N — > 00 pi 3 - One can then think of correc- 
tions to the continuum version of the height model for 
the CDM in terms of an 1 /N expansion (which we have 
not carried out). The results discussed here show that 
the 1/N corrections all the way down to N = 2 only 
correspond to a renormalization of the stiffness constant. 

VI. ORDER-PARAMETER DISTRIBUTION 

A columnar long-range ordered VBS on the square lat- 
tice breaks the translational and rotational lattice sym- 
metries. As we have seen in the previous sections, the 
RVB is a critical VBS with a rather slowly decaying 
dimer-dimer correlation function. This correlation func- 
tion, Eq. (jl.ip . measures the magnitude of the VBS order 
parameter. In this section we look at another aspect of 
these critical VBS correlations, probing the individual 
order parameters for columns forming with x and y ori- 
entation of the modulated bonds, defined as 

L L 

d x = z2(-^ x E^ y) ■ S (' T + i'^Uf, 
D v = E(- 1 )" E^ y) ■ s ^ y + 

y— 1 x—1 

where [...] CO nf indicates that these correlators are eval- 
uated for an individual configuration (i.e., in the RVB 
they are matrix elements between the sampled bra and 
ket states). The expectation values of these order param- 
eters vanish. In the CDM, the dimer-dimer correlation 
functions that we investigated before correspond to their 
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RVB W=(0,0) 



CDM W= (0,0) 



RVB All Sectors 



RVB W=(0,1) 



FIG. 16: (Color online) VBS order parameter distribution 
function P(D x ,D y ) in the space of point pairs (D x ,D y ) gen- 
erated in the Monte Carlo simulations of the RVB state (left) 
and CDM (right) for systems of size L = 64. We are only 
concerned here with the shapes of these distributions (a ring 
with depleted weight in the center for the RVB and a broad 
central peak for the CDM) and have therefore not labeled the 
graphs with the range of (D x , D y ) or the actual values of the 
probability density. 



squares, i.e., the dominant structure factor in reciprocal 
space (as seen in Fig. H]) is S(ir, 0) = (D%)/N, and the 
behavior of this quantity as a function of the system size 
is shown in Fig. GO In the RVB, as we have discussed in 
Sec. IIV Cl and Fig. [9j the squared order parameter based 
on the sampled values from Eq. (|6.ip is not exactly the 
same as the actual four-spin correlation function, but we 
have shown that the scaling properties are the same. 

We here study the probability distribution P(D x ,D y ) 
generated in the Monte Carlo sampling. Each generated 
configuration of the valence bonds corresponds to pair of 
values (D x ,D y ) evaluated according to the loop estima- 
tor Eq. (|2.12p . We use these to accumulate the histogram 
P(D X , D y ). Such histograms were generated by Suther- 
land in his loop-gas study^ and he noted a circular sym- 
metry of the distribution (instead of a 4-fold symmetry 
that would have been naively expected due to the lattice 
symmetry). At that time the results were affected by 
very large statistical uncertainties, however. 

Dimcr order-parameter histograms have recently be- 
come interesting in the context of deconfined quantum 
critical (DQC) point o 41 ' 42 in models exhibiting quantum 
phase transitions between the antiferromagnctic Neel 
state and a VBS state42di A Ion g-range ordered colum- 
nar VBS corresponds to a distribution P(D X , D y ) peaked 
at one of the four points (±|D|, ±|D|), with the magni- 
tude \D\ growing linearly with the system size N = L 2 . 
In a finite system, in which the Z4 symmetry is not bro- 
ken, one expects equal weight in all these four peaks, as 
well as some weight between the peaks (which is related 
to the tunneling probability between the four ordered 
VBS states). As a DQC point is approached from the 
VBS side, one expects an emergent U(l) symmetry in the 
system.— This is manifested in P(D X , D y ) as a circular- 
symmetric distribution , 42 ' 43 i.e., for a finite system size 
L, the discrete four-fold (Z4) symmetry naively expected 



FIG. 17: (Color online) VBS order parameter distribution 
function P(D X , D y ) for L = 48 systems in the grand-canonical 
winding number ensemble (left) and with winding number 
W = (0,1) (right). Here the range of D x ,D y values is the 
same in both cases, i.e., the distribution for W = (0, 1) is 
much narrower. 



for the VBS evolves into a continuous U(l) symmetric 
distribution. For fixed couplings, the Z4 symmetry devel- 
ops as L exceeds a length-scale characterizing the spinon 
confinement (which diverges at the DQC point). 

While the RVB is a critical state, it does not corre- 
spond to a DQC point, because the spin correlations 
decay exponentially. At a DQC point, both the spin 
and dimer correlations are critical4i It is nevertheless 
interesting to study the symmetry of the critical VBS 
order parameter in the RVB and to compare it with 
the corresponding distribution in the CDM [where in 
S(x,y) ■ S(x + l,y) is replaced by the dimer occupation 
number on the bond]. Results for L — 64 systems in the 
winding number sector W = (0, 0) are shown in Fig. 1161 
Completely circular-symmetric distributions are seen in 
both cases, with no signs of Z4 anisotropy. The natural 
expectation for a critical state is that the weight is cen- 
tered around (D x ,D y ) = (0,0), and this is in fact true 
for the CDM. Surprisingly, it is not true for the RVB crit- 
ical state: the distribution is instead ring shaped, with 
the dominant weight a finite radius away from the cen- 
ter. This is the behavior seen in candidate models for 
DQC points in the VBS state close to the phase transi- 
tion into the Neel state. The ring-shaped distribution in 
the RVB case is no contradiction to its being a critical 
state, because the ring's radius still grows slower with 
L than L 2 . The expectation value (D 2 )/N is twice the 
structure factor S(tt,0) and hence grows as L 2 ~ a , with 
a w 1.20 determined in Sec. IIV Cl 

In the case of a fixed non-zero winding number, the 
VBS order parameter is modulated by a plane wave, in 
the same way as its correlation function is, as discussed 
in Sec. IIV Cl Hence its spatial average tends to can- 
cel, with the result that the distribution function now 
has a central peak, as seen in Fig. [T7] (right panel) for 
W = (0, 1). For large winding numbers the distribution is 
marginally oval-shaped, reflecting the anisotropy induced 
by large winding numbers (see Appendix |B|) . In Fig. 1171 
the anisotropy is too small to observe clearly. Interest- 
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ingly, when all winding numbers are included in grand- 
canonical simulations, the ring-shaped distribution seen 
for W = (0, 0) in Fig.[16]no longer obtains. Although this 
sector completely dominates the grand-canonical ensem- 
ble (as seen in Fig. [5]), the narrow central peaks con- 
tributed by the non-zero winding number sectors com- 
pletely fill in the central portion inside the ring, resulting 
in a broad central peak, as shown in Fig. [TT] (left panel) . 



VII. SUMMARY AND DISCUSSION 

We have compared long-wavelength properties of 
short-bond RVB spin-liquid states with those of classi- 
cal dimers, specifically those associated with correlations 
and topological constraints of dimers. Taking properly 
into account the non-orthogonality of valence-bond basis 
states, arising from the internal bond-singlet spin struc- 
ture which is not present in classical dimers, we have car- 
ried out numerically exact Monte Carlo simulations of the 
four-point correlation function measuring the tendency 
to formation of a VBS state. In contrast to the exponen- 
tially decaying two-point spin correlations^ these VBS 
correlations decay as a power law. Such a power might 
have been anticipated based on the fact that the classi- 
cal dimcr-dimcr correlations decay as 1 /r 2 (although the 
overcompleteness of the RVB could in principle have led 
also to more dramatic deviations from the CDM), but the 
exact value of the exponent necessitates an exact treat- 
ment of the overcomplete basis, as we have done here. 
The result is that the correlations decay slower than what 
might have been anticipated, as l/r a with a « 1.20. 

The weighting of valence bond states is (qualitatively) 
different in that sampling the RVB state involves the 
transition graph of two states, whereas in the CDM only 
a single state is sampled (as different dimcr configura- 
tions are by definition orthogonal). In particular, the 
loops are small in the short-bond RVB, as they necessar- 
ily must be in order to give exponentially decaying spin 
correlations (whereas in an antiferromagnetically ordered 
state the typical loop size scales as the system siz o 17 ' 34 ). 
The operators that we measure are also different in the 
two systems: the "dimer-dimer" correlations in the RVB 
actually refer to two-spin operators, [Eq. (|1.2[) ] in place 
of just bond occupation numbers in the CDM. We have 
confirmed that the changed a exponent (and presumably 
other changed expectations) in the RVB state originate 
solely from the different state weighting, not from the 
form of the correlation- function estimator Eq. (|2.14[) . 

The RVB structure factor has a "pinch-point" at (tt, tt) 
in reciprocal space, in any winding number sector, like 
the well-known pinch-point in the CDM and other height 
models; it further shows singularities related to the crit- 
ical correlations near to (tt, 0) (but shifted by nonzero 
winding number) which are logarithmic for CDM at zero 
winding number, and otherwise are variable power laws. 
Finally, we found that introduced pairs of monomers, 
i.e., topological defects, are marginally (power- law) de- 



confined with a power law distribution of their separa- 
tions. 

Remarkably, all of the above observations fit into the 
framework of the "height model" with a stiffness constant 
K as worked out in Appendix |Bj Independent measure- 
ments of the stiffness constant can be derived from (i) 
logarithms of the probabilities of sectors with different 
winding numbers, (ii) the critical dimer correlation ex- 
ponent, (iii) the monomer pair separation exponent, and 
(iv) the pinch point of the structure factor 5(q). All 
yielded -FTrvb ~ 1-6-ffcDM- Other behaviors, which do 
not yield measurements of K, are also suggestive of this. 
Thus, our results vindicate at last the qualitative cor- 
rectness of the zero-overlap assumption adopted in the 
RK QDM, although quantitatively the RVB state has a 
larger degree of VBS order (as expressed by that ratio of 
stiffnesses 1.6). It is as if the RVB state were the ground 
state of the generalized RK state corresponding to some 
(still unknown) generalized classical dimcr model. 

We extended the model by introducing a small fraction 
of longer bonds (the next bipartite bond, which connects 
fourth- nearest neighbors). We studied the evolution of 
the power laws characterizing the dominant VBS corre- 
lations and monomer correlations as a function of the 
fugacity of long bonds. As in the CDM case, — , in the 
dimer-dimer correlations, a (tt, tt) modulated "dipolar" 
term continues to have the 1/r 2 behavior; on the other 
hand, a (tt, 0) modulated "critical" term has an increas- 
ing exponent, while the monomer- monomer distribution 
function has a decreasing exponent, both of which can 
be explained in terms of a decreasing stiffness for the 
"height" fluctuations. The monomers remain deconfined 
for all fugacities we studied. 

We further studied the modifications to correlations 
due to finite topological winding number, for both the 
RVB and classical dimers. The critical VBS correla- 
tions acquire a sinusoidal modulation, correlations be- 
come anisotropic, and the effective stiffness is increased, 
as expected from height-model calculations; 

We have also studied the joint probability distribution 
P(D x ,D y ) of the VBS order parameters for columnar 
order with x and y oriented bonds. We found this dis- 
tribution to be U(l) symmetric, which in analogy with 
the proposed deconfined quantum-critical points should 
correspond to the lattice-imposed Z4 symmetry of the 
VBS on the square lattice to be dangerously irrelevant 
[when regarded as a perturbation to an U(l) symmetric 
field theory] in these critical systems (both in the RVB 
and the CDM) . In a model that has one of these states as 
the ground state for some values of tunable parameters, 
e.g., the extended dimer models with "Cantor deconfine- 
ment" studied in Refs. [26| andHjJ one would then expect 
the U(l) symmetry to be emergent upon approach to the 
critical point. 

Although we have here studied the RVB state with- 
out reference to any specific Hamiltonian, some general 
conclusions can still be drawn based on our results. If 
a (local) Hamiltonian's ground state has algebraic corre- 
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lations, then it must correspondingly have gapless exci- 
tations. Thus, our results show that any Hamiltoniani^ 
with the RVB ground state is gapless in the singlet sec- 
tor, even though it has a spin gap. Furthermore, the close 
qualitative correspondence of the RVB static correlations 
to the RK model^ suggests the long-wavelength excita- 
tions are similar too; these are knowni to be coherent 
bosons with q dispersion. Some actual spin systems may 
be spin gapped but singlet gapless. This has long been 
claimed for the spin-1/2 kagome antiferromagnet ) 45 i 46 al- 
though the spin gap is small enough that an extrapolated 
value of zero can not be ruled out4£ From this viewpoint, 
it is interesting to verify that the original short-range 
RVB state has such a property. 

In experiments, the 2D organic 5 = 1/2 spin-liquid 
candidate, EtMe3Sb[Pd(dmit)2]2 shows gapless spin and 
singlet sectors in zero magnetic fieldj 4 ^ but in a magnetic 
field, spin excitations become gapped while singlet excita- 
tions remain gapless and have high mobility, as indicated 
by specific heat and thermal conductivity. 

On the theory side, one might ask whether our result 
should have been expected. Soon after the original pro- 
posal of the RVB wave function, field theorists argued 
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that it corresponded to a U(l) gauge theory 
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and 



for a "height model" to be in its rough phase, as we found, 
is equivalent to being asymptotically a (7(1) gauge the- 
ory. But, the numerical value of the stiffness constant 
K has not been measured previously (before our origi- 
nal estimate in Ref. [l8h : to our knowledge, it was not 
even suggested whether K should be larger or smaller 
than Kc dm of the QDM. If for no other reason, one 
must check the value of K since, were it much larger, one 
would find long-range order in the dimer correlations (a 
spin-Pcierls phase). 

It would clearly be interesting to try to derive the 
height model (or the continuum version of it) starting 
from an 1/N expansion of the classical dimer model, 
which corresponds to the RVB for SU(iV) spins in the 
limit N — > oo. Further, the recent construction^ of 
a model Hamiltonian which has exactly the RVB state 
studied here as its ground state also offers hope that 
one could actually, with extensions of that Hamiltonian, 
study a quantum phase transition in which the static 
properties of the critical point should be exactly those 
that we have investigated here in the RVB. 
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FIG. 18: (Color online) Action of a singlet projection operator 
in two different cases; (a) when the sites b, c are on different 
sublattices and (b) when b, d belong to the same sublattice. 
The arrows indicate the order of the spins in a singlet; (a, b) = 
(I ta-lb) — I 4-oti>))/v2, and, in the case of spins on different 
sublattices, conforms with the definition Eq. (|2.1[) of bipartite 
valence bond states. 



Appendix A: Four-spin correlators in the 
valence-bond basis 



In this appendix we work out the loop expression for 
four-spin correlators, analogous to the well-known two- 
spin expression Eq. (|2.12l) . 

It is useful to consider the singlet projectors 



Cij — — (Si • Sj — \). 



(Al) 



When acting on a valence bond, this operator is diagonal 
with eigenvalue 1. Denoting a singlet on sites a and b as 
(a, b), we have 



C ab (a,b) = (a,b), 



(A2) 



whereas acting on a pair of different valence bonds leads 
to a simple reconfiguration of those bonds, e.g., 

C bc (a,b){c,d) = i(c,b)(a,d), (A3) 
C bd (a,b)(c,d) = ±(a,c){b,d), (A4) 

which can be shown easily by going back to the basis of f 
and i spins. Note the order of the indices within the sin- 
glets in Eq. (|A3[) . which reflects consistently the chosen 
convention in the valence-bond state definition Eq. f|2 . 1 [) 
when the sites a, c are on sublattice A and b, d on sublat- 
tice B. We will also have to consider operations on two 
spins belonging to the same sublattice, as in Eq. (|A4|) . 
We have not specified a convention for the order of the 
spins in singlets formed between two spins on the same 
sublattice, therefore, it is important to keep track of the 
signs, which depends on the order in which the singlets 
are written. 

Figure [T5] illustrates the two different types of sin- 
glet projector outcomes in Eq. (|A3|) and Eq. (|A4I) . In 
Fig. [TSta), both the initial and the final bond pairs are 
bipartite whereas in Fig. IT8T b) the bonds after the op- 
erator has acted are non-bipartite. The non-bipartite 
bonds do not belong to the restricted basis of bipartite 
valence-bond basis in which we normally work. However, 
when generating non-bipartite bonds such as this (which 
can happen in the course of calculations), we can always 
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FIG. 19: (Color online) Illustration of the equivalence 
Eq. ()A5|) , due to overcompleteness, between a state formed 
by two non-bipartite valence bonds and a superposition of 
two states involving only bipartite bonds. 

rewrite them in terms of bipartite bonds. One can eas- 
ily verify the following equivalence between valence bond 
pairs; 

(a, c) (b, d) = (a, b) (c, d) - (a, d) (c, &) , ( A5) 

which is illustrated in Fig. [19] This relationship is par- 
ticularly useful when sites a, c € A and b,d € -B, but it 
of course holds irrespective of sublattices. 

As in Eq. (|2.11[) . we can take advantage of the spin- 
rotational symmetry also when considering a four-spin 
correlation function, writing the corresponding matrix el- 
ement as 

(V f3 \(S k -S l )(S i -S j )\V a )=3(Vp\S* k S?(S i -S j )\V a ). (A6) 

Note, however, that we cannot further reduce this ex- 
pression to a correlation function involving only z-spin 
components, because if 7 7^ z, 

(Vp\S z k S?S?S*\V a ) ± (Vp\S* k SfS2S]\V a ). (A7) 

It is easy to see that the matrix element Eq. (|A6|) is non- 
zero only if all four indices i,j,k,l belong to the same 
loop, or if there are two indices in each of two loops. To 
carry out the calculations for these cases, it is convenient 
to make use of the singlet projection operator Eq. (|A1[) 
and write the matrix element as 

(V fS \S* k Sf(S i -S j )\V a }= (A8) 
\(V p \S* k S?\V a ) - (V \S z k S?dj\V a ). 

We only go through the calculation in detail for the case 
where all four indices belong to the same loop, which is 
the most complicated situation. 

The procedure is illustrated in Fig. [20] Acting first 
with the singlet projector dj, the loop is split into two 
separate loops if i,j are on different sublattices, as shown 
in Figs. I2"07a) and l2"07 b). If these sites are on the same 
sublattice, as in Fig. [207 c). the loop instead becomes 
"twisted" by two non-bipartite bonds. This loop can be 
re-cast in terms of two different contributions containing 
only bipartite bonds, by using the valence-bond equality 
illustrated in Fig. QjU In each case, after dj has acted, 
we can return to the spin representation of the valence 
bonds and evaluate the average of the remaining opera- 
tor S k Si exactly as we did for the two-spin correlation 
function. Here the result depends on whether k, I are in 
the same loop (giving a non-zero correlation) or differ- 
ent loops (giving a zero average) after the loop-splitting 
with dj has been enacted; these two different cases are 




FIG. 20: (Color online) Operations for evaluating the four- 
spin matrix element (V8|(Sfe-S;)(Si-Sj)jV Q ) when all the sites 
i,j, k, I are in the same loop of the transition graph. The thin 
lines connecting labeled sites refer to the operator components 
SfcSf and dj in Eq. (|A8[) . The solid and dashed bonds belong 
to \V a ) and {Vp\, respectively. 

illustrated in Fig. [2DTa) and [207 b) for the case i, j in dif- 
ferent sublattices [while for i, j on the same sublattice, 
Fig. I2"07c) only shows the case of k, I in different parts 
of the split loop]. In all cases, the matrix element ratio 
(Vp\S k Sidj\V a ) /(Vp\V a ) is now easy to compute using 
Fig. [50] and keeping in mind that an increased number 
of loops after a split by CV,- increases the corresponding 
matrix element by a factor 2 according to the loop expres- 
sion Eq. (|2.6p for the overlap. The four-spin correlation 
can then be extracted using Eqs. (|A6[) and (|A8[) . 

In order to write the final result in a compact unified 
form for all the different cases, it is useful to introduce 
the concept of subloops with respect to the operator CV,- 
of a loop containing sites or (i, j)-subloops. As seen 
in Fig. [20J regardless of whether i,j are on the same or 
different sublattices, the loop is split in the same way by 
Cij in all cases where such split loops appear. This can be 
formalized by the following convention: The splitting of a 
loop into (i, j)-subloops is accomplished using the bonds 
in the ket \V a ) (the solid bonds in Fig. [20] on which dj 
acts), i.e., the two V^-bonds on which i,j are located are 
those that are reconfigured in such a way that the loop 
splits into two. The subloops then always contain only 
bipartite bonds. This definition is illustrated in Fig. [21] 
We also introduce a symbol to distinguish between the 
cases of k, I in the same subloop or different subloops; 

$ k ! — J ^' ^ or ^' ' m ^ e same j)-subloop, (A9) 
13 [ 1, for k, I in different (i, j)-subloops. 

If i,j are on the same bond of \V a ), dj does not change 
the loop and there is then only a single subloop (the 
intact original loop) and 6fj = for all k, I. 

The remaining cases of non-zero four-spin matrix ele- 
ments involve two loops (with two indices in each loop). 
These calculations are easier than the case of all indices 
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FIG. 21: (Color online) Subloops of a valence-bond loop with 
respect to two sites The "cuts" splitting the loop into 
subloops are at the solid bonds connected to i and j (which 
belong to the ket \ V a ); the state on which dj acts), irrespec- 
tive of the two possible locations of i,j within these bonds. 
When i,j are sites in the same bond in \V a ), there is only a 
single subloop (the whole loop). 

in the same loop, because there are no subloops to con- 
sider, and we just list the results. The full final result 
for all non-zero four-spin matrix elements is given in the 
main text as Eq. (|2.14p . 

Note that whereas the sign of the two-spin correlation 
Eq. (|2.12|) is always dictated by the staggered phase fac- 
tor, the sign of the four-spin correlation is different from 
the four-site staggered phase 4>ij4>ki if all the indices are in 
the same loop and k,l belong to different (i, j)-subloops. 

The concept of subloops may seem unnecessarily com- 
plicated in the definition of 6% in Eq. (|A9[) . since this 
number (0 or 1) is essentially also determined by the or- 
der in which the indices i,j,k,l appear when traversing 
a loop. If only one of the indices fc, / appear between 
then, in most cases, k,l are in different subloops 
and 5 % i, = 1. There are, however, special cases where 
the definition based on the order of indices is ambiguous, 
e.g., when they are all on the same valence bond in the 
ket \V a ). In that case, fc, I are in the same subloop and 
<5£J = 0, as also explained in Fig. [3T] 

Appendix B: Calculations based on height 
representation 

Any complete covering of a bipartite planar lattice 
(such as the square lattice) by dimers can be mapped 
into a configuration of "heights" representing a kind of 
interface model. Often , the ensemble weighting corre- 
sponds to the "rough" phase of the interface. In this 
case, many statistical properties may be derived from a 
simple (Gaussian) classical field theory in terms of the 
coarse-grained height function, using the "Coulomb-gas" 
formalisms introduced in the Kostcrlitz-Thoulcss theory 
of the two-dimensional XY modeh 49 ' 50 Bipartite dimcr 
coverings are a subset of a larger class of "height" mod- 
els treated by this formalism, which also include random- 
tiling quasicrystalsi 39 ' 40 

The CDM is known to be in this "rough" phase. In the 
case of the RVB wave function, for which this property 
had not been known, it is shown in this paper that all 
statistical behaviors are consistent with a rough height 
model. It should be emphasized that this is an emer- 
gent behavior, since there is no exact way to map spin 



states to dimer coverings (the dimers to spins mapping 
is not invertible). We might hypothesize the existence 
of some hidden, nonlocal way to define winding num- 
bers and perhaps height fields from the spins; however, 
the nonzero overlap between configurations in different 
winding- number sectors (see Fig. [5]) shows that there can 
not be an exact mapping of that sort. 

The starting point of the height treatment is that the 
probability of a (coarse-grained) height field {h(r)} is 
given by exp[-F tot ({h(r)})], where 

F tot = J d 2 r±K\Vh(r)\ 2 . (Bl) 

We here study various consequences following from this. 

1. Relation of height field and dimer operators 

There are two closely related ways to define a height 
function, for a dimer model, as laid out in Ref. H3- The 
microscopic height h(r) is defined on dual vertices (cen- 
ters of plaquettes); we set 

h{x+\,y+\) - h(x-\,y+\) 

= (-ir + y[4n y (x,y)-l], (B2a) 
h(x + \,y+\) - h{x + \,y-\) 

= (-l) x+ y[4n x (x,y)-l]. (B2b) 

Thus h takes a step ±3 across a dimer, or =pl across 
an unoccupied bond, where the sign alternates between 
even and odd vertices of the lattice. If one takes four 
steps around a vertex, one crosses a dimer once and an 
unoccupied bond three times such that the net difference 
is zero, ensuring a well-defined height field. 

A second, locally averaged height function h(x, y) is de- 
fined on the original vertices, being the mean of h on the 
four surrounding plaquettes. [Note the locally averaged 
h(x,y) is not quite identical to the fully coarse-grained 
height function assumed in the field theory, although we 
use the same notation h(r).] This h(x,y) is uniform in 
any one of the four special domains in which the dimers 
are aligned on opposite sites of plaquettes; it shifts by 
one unit on crossing a domain wall to the next domain. 
A change of ±4 in h brings us back to the same domain. 

Thus, the dimer occupation can be written as a period- 
four function of the local height variable, 

n x (r) = i [cos (^) 2 + (-1)* cos (^)], (B3a) 

% (r) = I[sin(^) 2 + M) y sin(^)]. (B3b) 

The configurations with a given winding number may be 
visualized as fluctuating domains with smoothed domain 
walls. For winding number W = (W x ,0), a net num- 
ber of domain walls 4W X must be crossed as the system 
is traversed in the x direction. There is no long-range 
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dimer order, so the domain walls thereby enforced are 
delocalized; indeed, in a snapshot of the configuration, 
they are lost in the dense array of random domain walls 
which arc part of the inherent fluctuations even in the 
W = (0, 0) sector. 



2. Effects of long dimers 

In the present simulations, sometimes dimers are per- 
mitted (both in CDM and RVB models) between sites 
separated by a (2, 1) type vector with a fugacity Z 2 . This 
requires us to modify the height construction. Say this 
dimer extends from (0,0) to (2,1). The height changes 
across the lattice edges (0,0)-(l,0) and (1,1)— (2,1) as if 
there were ordinary dimers occupying both edges (i.e. 
— 1 times the height change if those edges were vacant.) 
As for the lattice edge (1,0)— (1,1) bisected by the long 
dimer, the height change is +5 times the height change 
the vacant edge would have had. Around the vertex (1,0) 
or (1,1), the net height changes are 3 + 3 — 5 — 1 = 0, 
showing the modified construction is well defined. 

It can be seen that long dimers allow larger differences 
in height between adjacent sites. In the coarse-grained 
picture this means that height gradients are penalized 
less and thus K is decreased. Indeed, it was observed in 
previous worki^ that in the CDM when only long dimers 
are present, K is reduced by a factor of 2/9. 



3. Dimer correlations: dipolar term 



gradient terms in Eq. (|B4[) . we find 



It seems as if Eqs. (|B2j) and (|B3|) express contradictory 
relations between the height field and the dimer configu- 
ration. The proper resolution is that the dimer field has 
two slowly varying parts that are modulated in different 
ways with respect to the lattice, 
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which is equivalent to Eq. (2.4) of Ref. [2fJ. It turns out 
that the n x -n x dimer occupation correlation, as a func- 
tion of displacement r = (x, y), breaks up into two slowly 
decaying terms, D xx (r) = _D^ p (r) + _D™'(r), which are 
due to the two kinds of terms in Eqs. (|B4|) . 

Consider the first kind of term. Equation (|B1[) implies, 
for the Fourier transform of the height field, (|ft,(q)| 2 ) « 
1/A'|q| 2 for small wavevectors q. Combining with the h 
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for the x-dimer structure factor near Q = (%, tt). Tak- 
ing the Fourier transform of Eq. (|B5|) gives the (two- 
dimensional) pseudo-dipolar correlations 

2 2 

D^ir) « (-l) a+ »Const ^ ~ V U . (B6) 
xx v i \ i 27rA|r| 4 v ; 

The radial dependence of this is 1/r 2 in any direction, 
irrespective of the value of K. 



4. Dimer correlations: Critical term 

We now turn to the second kind of term in Eqs. (|B4I) . 
the terms periodic in h. By a calculation standard in 
height-model literatu re 37 ' 38 , they imply the Coulomb gas 
(critical) term, 



D x f(r) cx 



(-1)" 



where 



a 



2ttA 



7T 

8K' 



(B7) 



(B8) 



It is a peculiarity of the CDM, with nearest-neighbor 
dimers and equally weighted configurations, that a = 
2. Thus both terms have the same decay exponent and 
in fact they cancel exactly on certain sites. Modifying 
the relative weighting of dimer configurations normally 
changes a. If a < 1/4, the height configuration locks into 
a flat state (roughening transition) which means that the 
dimers lock into a long-range ordered state. However, in 
this study, a is reduced from the CDM value of 2 by a 
relatively modest amount. 

The same kind of calculation implies that 



1 



flS'(L/U/2) K - 



(B9) 



with the same a as in Eq. (|BT|) . but a different pref- 
actor. Note that (so long as the elasticity is isotropic) 
the dipolar contribution Z?^ p (r) is exactly zero along 
the lines x = ±y (even as its asymptotic r depen- 
dence breaks down) and therefore does not contribute 
to D xx {L/2,L/2). 



Topological (monomer) defects and their 
correlations 



If a site is uncovered, the height differences do not can- 
cel in going around it, but change by b — ±4 (where the 
sign depends on whether the vertex is even or odd). Such 
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defects can only be created in pairs of opposite charge, 
and play the same role as vortices in the Kostcrlitz- 
Thoulcss theory. The K values in our simulations are 
small enough that we are above the Kosterlitz-Thouless 
unbinding transition, i.e., if there were nonzero fugacity 
to have defects, they would destroy the critical state at 
sufficiently long length scales. However, the fugacity is 
in fact zero (except that in some simulations, one pair is 
inserted by hand as a probe). 

The presence of a defect at (say) the origin enforces 
a background gradient in the height field with \Vh\ = 
b/2irr. When substituted into Eq. (|B1|) . that would give 
a logarithmically divergent total, except that the diver- 
gence gets cut off by another defect at distance R. The 
result is that the effective potential cost for the defects 
to be separated by R is (K/2ir)b 2 InR, and the pair dis- 
tribution is given by 



with 



Kb 2 8K 



(BIO) 



(Bll) 



2-7T 7T 

and in particular /? = 1/2 for the basic CDM. 

6. Sector probabilities 

We now turn to the effects of enforcing net wind- 
ing numbers W x ,W y . This is equivalent to a boundary 
condition that h(L,y) = h(0,y) + AW X and h(x,L) = 
h(x, 0)+4W y . In light of Eq. (|B3[) . no discontinuity is im- 
plied in the actual dimer pattern, since that depends on 
h(r) with period 4. It would be exactly analogous to en- 
forcing, in an XY model, angle differences (2irW x , 2irW y ) 
across the system. 

Thus the effect of winding number (W x , W v ) is to im- 
pose a uniform "background" height tilt (m x ,m y ) = 
i(W x ,W y )/L. We write 



h(r) = m x x + m y y + h'(r), ) 



(B12) 



separating the height field into the background plus a 
(smaller) deviation h'(r) that satisfies periodic boundary 
conditions. 

If we substitute the free energy Eq. (|B1[) into 
Eq. (|B12[> . we see that 

F tot (fh}) = F tot {{h'}) + AF(W X , W y ), (B13) 

where 

AF(W x ,W y ) = \KL 2 (ml + w 2 ) = 8K(W 2 + W 2 ). 

(B14) 

Since F tot in Eq. (|B13[) is exactly the same function 
as before, it follows that when we integrate over all 
configurations of {h'(r)} to obtain the partial partition 



function Z(W x ,W y ) for a given sector, Z(W x ,W y ) — 
Z(O i O)exp[-AF(W x ,W y )}. We conclude that the rel- 



ative probabilities of different sectors are given by 



P(W x ,W y ) = P(0,0)e 



(B15) 



In checking the normalization of P(W x ,W y ), it should 
be remembered that e.g. the (1,0) sector is fourfold de- 
generate [the possible winding numbers are (±1,0) and 
(0,±1)], as are the (1,1) and (2,0) sectors. 



7. Correlation modulation due to winding number 

To calculate the critical contribution in the presence 
of a background h gradient associated with a winding 
number, we merely need to substitute Eq. (|B12[) into 
Eqs. (|F34|) . remembering that the rightmost terms are the 
ones contributing to the desired correlation. The result 
is that we get the correlation due to the h' field (i.e. the 
same as before) times cos[^(m x a: + m y y)], where (x,y) 
is the vector connecting the two points. In other words, 



^ r i t (r;^)=I?« it (r;0)cos( ( 5Q-r), 



(B16) 



where D™*(r; W) means D™*(r) given winding numbers 
W, and 



SQ 



2tt 



(m x ,m v ) = 2n(W x ,Wy)/L. 



(B17) 



Since £)£"*(!■; 0) already includes a (— l) x modulation, it 
follows that the structure factor singularity of D" a ! t (r; W) 
gets shifted to 



Q = (tt,0)±(5Q. 



(B18) 



8. Anisotropic effects due to winding number 

In a height model, the free-energy density is a function 
of Vh(r) and its derivatives, satisfying all lattice symme- 
tries. The free-energy density in Eq. (|B1|) is the lowest 
term of its Taylor expansion in \7h. The next terms 
consistent with the square lattice are quartic, thus, the 
free-nergy density becomes 



f(Vh) = ±K\Vh\ 2 + 911 



\dx/ \dy ) 



( dh\ 2 ( dh\ 2 



If we insert Eq. (|B19| into Eq. (|B12|) The effective free 
energy density to lowest order in b! is 

1 (dh\ 2 1 (dh\ 2 (dh\ (dh\ ,„ . 

f=M — ) +-K y {—) +K X y( — ){—), (B20) 



' \dx 



\dy 



V dx I \dy 
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where 



A", 

Ky 



X + 12g n ml + 2g 12 ml 



K + 12g xl m\ + 2g 12 i 
4g 12 m x m y . 



(B21a) 
(B21b) 
(B21c) 



The nonlinear terms of a background tilt were consid- 
ered and measured from simulations in the quasicrystal 
random tiling context^. It is possible, in principle, to ex- 
tract analytical expressions for the nonlinear terms from 
the exact solutions. 

Next we consider how this modifies correlations. For 
simplicity, consider the case m y — 0. We make a change 
of variables 



x' = 7x; y' = 7 1 y, 



(B22) 



where 



7 = {K x /K y f/±. (B23) 
In the new coordinates, the free energy density is 



/ = ~ 2 K< 



\dx i ) + \~dy~ 7 



(B24) 



with an effective stiffness K' = yJK x K y , In these new 
coordinates, Eq. (|B24|) looks isotropic again and the same 
results must follow for the behavior of all correlations. 
In particular, the dimer and monomer correlation decay 
exponents, a and /?, depend on K' in the same way they 
previously did on K. In the general case that m x m y =/= 0, 
the effective stiffness is 



K' 



K X K, 



y 



K 2 

xy 



(B25) 



For small W/L, i.e. small (m x ,m y ), this reduces in light 
of Eqs. (JB2U to to K> m K + 96( gil + g 12 )(W 2 + W 2 )/L 2 . 
Hence large L, and a winding number W the corrections 
to exponents scale the same way, 5a ~ 5(3 ~ W 2 /L 2 . 

Notice that the decay exponent is the same in all spa- 
tial directions. The way the anisotropy gets expressed 
in the correlations with variable exponents is that (e.g.) 
dimer correlations do not fall off exactly as l/r a , but 
as 1/r , 



V7^ 



■r 



~ 2 y 2 , and simi- 



rather as 1/r' , where r 
larly for monomer pair separations. It would be interest- 
ing to see whether the anisotropy of spin correlations, as 
shown in Fig. [51 is expressed by the same ratio 7. 
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